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Ly-a RADIATIVE TRANSFER IN COSMOLOGICAL SIMULATIONS AND APPLICATION TO A 
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ABSTRACT 

We develop a Ly-a radiative transfer (RT) Monte Carlo code for cosmological simulations. High res- 
olution, along with appropriately treated cooling, can result in simulated environments with very high 
optical depths. Thus, solving the Ly-a RT problem in cosmological simulations can take an unrealisti- 
cally long time. For this reason, we develop methods to speed up the Ly-a RT. With these accelerating 
methods, along with the parallelization of the code, we make the problem of Ly-a RT in the complex 
environments of cosmological simulations tractable. We test the RT code against simple Ly-a emitter 
models, and then we apply it to the brightest Ly-a emitter of a gasdynamics-|-N-body Adaptive Refine- 
ment Tree (ART) simulation at z ~ 8. We find that recombination rather than cooling radiation Ly-a 
photons is the dominant contribution to the intrinsic Ly-a luminosity of the emitter, which is ~ 4.8 x 10*^ 
ergs/s. The size of the emitter is pretty small, making it unresolved for currently available instruments. 
Its spectrum before adding the Ly-a Gunn-Peterson absorption (GP) resembles that of static media, 
despite some net inward radial peculiar motion. This is because for such high optical depths as those 
in ART simulations, velocities of order some hundreds km/s are not important. We add the GP in two 
ways. First we assume no damping wing, corresponding to the situation where the emitter lies within the 
HII region of a very bright quasar, and second we allow for the damping wing. Including the damping 
wing leads to a maximum line brightness suppression by roughly a factor of ~ 62. The line fluxes, even 
though quite faint for current ground-based telescopes, should be within reach for JWST. 

Subject headings: cosmology: theory — diffuse radiation — galaxies: formation — intergalactic 

medium — radiative transfer — line: formation — radiative transfer — resonant — 
polarization 



1. INTRODUCTION 

Since the classic paper by Partridge & Peebles (1967), 
intense observational efforts have focused on the search 
for Ly-a emitters at high redshifts. Although most of 
the early attempts ended in negative results before the 
mid 1990s, recent observational advances enabled us to 
identify star forming galaxies at ever increasing redshifts. 
Currently, several observational projects, such as LALA 
(e.g., Rhoads et al. 2003), CADIS (e.g., Maier 2002), the 
Subaru Deep Field Project (e.g., Taniguchi et al. 2005), 
etc., spectroscopic surveys that use lensing magnification 
from clusters (e.g., Santos et al. 2003), surveys that com- 
bine Subaru (e.g., Hu et al. 2004) or HST/ACS/NICMOS 
imaging (e.g., Stanway et al. 2004; Dickinson et al. 2004, 
for the GOODS survey) with Keck spectroscopy, etc., focus 
on finding high-z starforming galaxies. Surveys currently 
reach up to z ~ 7 — 8 (e.g., see Bouwens et al. 2004, for re- 
cent results from the NICMOS observations of the HUDF) 
and will likely reach higher redshifts in the coming years 
(e.g., via JWST). 

The hydrogen Ly-a line is a very promising way to probe 
the high-redshift universe. Besides yielding redshifts, the 
shape, equivalent width and offset of the Ly-a line from 
other emission/absorption lines potentially convey valu- 
able information about the geometry, kinematics, and un- 
derlying stellar population of the host galaxy. Further- 
more, after escaping the environment of the host galaxy, 
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Ly-a photons are scattered in the surrounding IGM. The 
presence or absence of observed Ly-a emission can be used 
to place constraints on the state of the IGM, useful in con- 
straining for example the epoch and topology of reioniza- 
tion. Because of the numerous factors that contribute to 
the final Ly-a emission, the interpretation of such features 
can be very complex. To use all the currently available 
and future observations in the most effective way possi- 
ble we need to improve our theoretical understanding of 
Ly-a emission from high-redshift objects. To this end we 
develop a general Ly-a radiative transfer (RT) scheme for 
cosmological simulations. As an example, we apply the 
RT scheme to gasdynamics-|-N-body Adaptive Refinement 
Tree (ART; Kravtsov 2003) simulations of galaxy forma- 
tion. 

There are quite a few studies of Ly-a emission from 
high-z objects (e.g., Gould & Weinberg 1996; Haiman & 
Spaans 1999; Loeb & Rybicki 1999; Ahn et al. 2001, 2002; 
Zheng & Miralda-Escude 2002; Santos 2004; Dijkstra et al. 
2005a, b). The problem these studies address is highly 
complex and has many unknowns. Inevitably, most of 
these studies had to make at some point some simplify- 
ing assumptions. Usually, a high degree of symmetry for 
the emitting source, and its density, temperature, and ve- 
locity field is assumed. The same is the case with respect 
to the processes that are responsible for the production of 
Ly-a photons. On the other hand, cosmological simula- 
tions hopefully capture most of the basic elements, lifting 
thus practical constraints that existed in these previous 
studies. 

There is a small number of related studies using cos- 
mological simulations (Fardal et al. 2001; Furlanetto et al. 
2003; Barton et al. 2004; Furlanetto et al. 2005; Cantalupo 
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et al. 2005; DcUiou ct al. 2005a, b). Some of these simula- 
tions are lacking crucial processes such as radiative cooling 
of the gas and consistent RT, the various sources of Ly-a 
photons, and/or sufficient resolution in order to resolve the 
dumpiness of the gas. Furthermore, most of these studies 
do not perform Ly-a RT, but rather they assume that the 
observer sees whatever is being emitted initially, simply 
modified by e~'^ with r the optical depth for Ly-a scat- 
tering due to neutral hydrogen between the emission point 
and the observer. Namely, in most cases Ly-a spectra from 
simulations are treated as absorption spectra when, in real- 
ity, they are scattering spectra (see, e.g., Gnedin & Prada 
2004). For gas well outside the source of emission this 
is an appropriate approximation since scattering off the 
direction of viewing removes the photons that could be 
observed and thus appears as effective absorption. This 
is no longer true for the source of emission itself, since 
photons that were originally emitted in directions differ- 
ent from the direction of observation may scatter into this 
direction. 

It is important that the difficulty of implementing a Ly- 
a RT scheme for cosmological simulations become clear. 
The classical problem of resonance RT, relevant to a wide 
range of applications from planetary atmospheres to ac- 
cretion disks, has been quite extensively studied in the lit- 
erature (e.g., Zanstra 1949; Hummer 1962; Auer 1968; Av- 
ery & House 1968; Adams 1972; Harrington 1973; Neufeld 
1990; Ahn et al. 2001, 2002). However, analytical solu- 
tions derived in the past are applicable only to certain 
specific conditions. On the other hand, the slow conver- 
gence of the numerical techniques used limited the nu- 
merical studies at optical thicknesses that are relatively 
low compared to those encountered in high-rcdshift galax- 
ies (and cosmological simulations of high-redshift galaxies, 
as we will show). Thus, unlike previous studies, most of 
which focused on the classical problem of resonant RT in 
a semi-infinite slab, in cosmological simulations one has to 
solve simultaneously thousands or even millions of these 
problems. ^Furthermore, having in mind existing and fu- 
ture cosmological simulations that can achieve sufficiently 
high resolution to resolve the gas dumpiness and that treat 
cooling appropriately, we anticipate column densities that 
are orders of magnitude higher than those found in lower 
resolution simulations without cooling. In this case we 
need a RT algorithm much faster than the more standard 
direct Monte Carlo approach [which, however, is our start- 
ing point] of previous studies. Thus, we must develop RT 
acceleration methods that, along with the highly parallel 
nature of thc> RT problem that enables us to make use of 
many parallel machines, can make the Ly-a RT problem 
tractable. 

The paper is organized as follows. In § 2 we discuss the 
RT scheme. More specifically, in § 2.1 we present the ba- 
sic Monte Carlo algorithm, in § 2.2 we present tests of the 
basic algorithm, in § 2.3 we discuss the acceleration meth- 
ods we use to speed up the RT, and in § 2.4 we present 
the method images and spectra are constructed. In § 3 
we discuss in detail an application of the Ly-a RT code to 
ART simulations. More specifically, in § 3.1 we briefly give 

^ For example, in the case of Adaptive Mesh Refinement (AMR) 
codes, each time a photon enters a simulation cell one has the equiv- 
alent of a new slab RT problem. 



some information about the ART simulations. In § 3.2 we 
discuss the intrinsic Ly-a emission of the specific simu- 
lated Ly-a emitter we focus on. In § 3.3 we present results 
on the emitter before RT. In § 3.4 we discuss results after 
performing the RT, and with/ without the Gunn-Peterson 
(GP) absorption, as well as with/without the red damp- 
ing wing of the GP absorption. In § 4 we discuss and 
summarize our results and conclusions. 

2. THE LY-a RADIATIVE TRANSFER 
2.1. The basic Monte Carlo code 

The following discussion assumes in various places a cell 
structure for the simulation outputs, as is inherently the 
case in AMR codes. However, the Ly-a RT code we dis- 
cuss is applicable to outputs from all kinds of cosmological 
simulation codes, since one can always create an effective 
mesh by interpolating the values of the various physical 
parameters. The size of the mesh cell can be motivated by 
resolution related scales (e.g., the softening scale, or larger 
if convergence tests with respect to the Ly-a RT justify a 
larger scale). Thus, in what follows we refer to simulation 
cells either the direct output of the cosmological simula- 
tion code has a cell structure or not. 

The initial emission characteristics (simulation cell, fre- 
quency, etc.) of each photon depend on the specific phys- 
ical conditions, thus we defer this discussion for §3 where 
an application to a Ly-a emitter produced in ART cos- 
mological simulations is presented. After determining the 
initial characteristics for each photon, we follow a series 
of scatterings up to a certain scale where the detailed RT 
stops. This scale is to be determined via a convergence 
study. In this subsection we describe the basic steps of 
the algorithm. 

2.1.1. Propagating the photon 

For every scattering we generate the optical depth, which 
determines the spatial displacement of the photon, by sam- 
pling the probability distribution function e~'^ 

T = -\n{R), (1) 
with R a uniformly distributed random number. This op- 
tical depth is equal to 
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with uhi the number density of neutral hydrogen. The 
function ctl is the scattering cross section of Ly-a photons 
as a function of frequency, defined in the rest frame of the 
hydrogen atom as 

=2 Ai/i/27r 
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where /12 = 0.4162 is the Ly-a oscillator strength, z/q = 
2.466 X 10^'"' Hz is the line center frequency, Ai^l = 9.936 x 
10^ Hz is the natural width of the line, and other symbols 
have their usual meaning. In equation (2) the fact that the 
photons are encountering atoms with a Maxwellian distri- 
bution of thermal velocities has been taken into account. 
Integrating over the distribution of velocities, the resulting 
cross section in the observer's frame is 



o-(a;) = /12 
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where 



H{a,x) = - r 



{x — y)"^ + o? 



dy 
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is the Voigt function, x = {v — vq)I l\vu is the relative 
frequency of the incident photon in the observer's frame 
with Ai^£) = yJ^ksT I {m,pc'^)iyo the Doppler width, and 
a = Ah'L/2Av£i with Ai^l the natural line width. Assum- 
ing that a is independent of /, the optical depth is given 

by 

T = nHi(j{x)l . (6) 

When applied to cosmological simulations, equation (6) is 
substituted by a sum of terms similar to the r.h.s. This 

sum is over the different cells (=different physical condi- 
tions such as neutral hydrogen density, temperature, etc.) 
that the photon crosses until it reaches r and gets scat- 
tered. 

For the Voigt function we use the following analytic fit, 

which is a good approximation to better than 1% for tem- 
peratures T > 2K (N. Gncdin, personal communication) 

= -lAi H{a,x) = -r^^(a;) 
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where i = x^, and g = if z = (i - 0.855)/(i 3.42) < 
and 



q = z{\ 
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7r(x + 1) 

X {0.1117 + z [4.421 + z(-9.207 + 5.674z)]} (8) 

if 2; > 0. The definition in terms of the function 4>{x) is 
also given since the latter has been used in many previous 
studies, and we also use it in what follows. If in addition to 
the thermal motion of the atoms there is bulk motion, such 
as peculiar or Hubble flow velocities, in equation (4) we use 
Xf = X — (^'/z/c)^'o/A^'£), where Vfz is the component of 
the fluid bulk velocity along the direction of the incident 
photon. 

In equation (2) the cross section a becomes /-dependent 

when Hubble expansion is taken into account. In this case 
the equation is an integral and does not reduce to the sim- 
ple algebraic equation (6). To propagate the photon one 
must solve for the step which is the upper limit of the in- 
tegral. In the simple examples discussed in §2.2, things 
arc relatively simple even when the Hubble expansion is 
included, since in these cases there is homogeneity and 
isothermality and no sum over cells is required. In those 
cases, Hubble expansion is included as follows: 1.) we 
make a first guess for I using the Hubble velocity at the 
current point, 2.) we use as a step for the photon a certain 
fraction of 3.) for a specified tolerance with which we 
want to achieve r, we refine the step as necessary. Note 
that the simple tests of the code presented in §2.2 do not 
include peculiar motions. In the actual simulations the 
peculiar velocities rather than the Hubble expansion are 
dominant on the relevant scales (e.g., for the emitter we 
focus on, the mean radial component of the peculiar mo- 
tion dominates over the Hubble expansion up to about 80 
physical kpc). In the detailed RT which we perform within 
such distances, we approximate the subdominant Hubble 



expansion velocity within a certain cell by the expansion 
velocity that corresponds to the center of that cell. This 
is calculated to have a negligible effect on the results. 

The n = 2 state of atomic hydrogen consists of the 
251/2; 2P1/2 Etnd 2P3/2 substates, whereas the n = l state 
consists of 15'i/2- According to the electric dipole selec- 
tion rules, the allowed transitions are 2P1/2 to 15i/2 and 
2P3/2 to 151/2, whereas 25i/2 corresponds to destruction 
of the initially absorbed Ly-a photon, since this state de- 
excites throiigh the emission of two continuum photons. 
The multiplicity of each of these states is 2J -|- 1. Thus 
the probabilities for the 2P states the atom can be found 
in when absorbing the Ly-a photons 2Pi/2 : 2P3/2 are 
1:2. Collisions can potentially cause the 2P — s- 2S tran- 
sition in which case the photon gets destroyed. A similar 
destruction effect can be caused by the existence of dust. 
Both these destruction mechanisms are briefly discussed in 
the context of the ART simulations in §3.5.1 and §3.5.2, 
respectively. Considering the 2Pi/2 and 2P3/2 cases sepa- 
rately, one would have to modify both the Voigt function 
and the velocity distribution of the scattering atom dis- 
cussed in the next section (see, e.g., Ahn et al. 2001). 
However, the level splitting between the two 2P states is 
small, just 10 GHz. This corresponds to a velocity width 
of ^ 1 km/s, much smaller than the width due to ther- 
mal velocities in media with roughly T > 100 K. In addi- 
tion, even for lower temperatures, this level splitting is still 
small for high optical depths. In our case, the thermal, pe- 
culiar, and Hubble velocities are all more important than 
the splitting, and combined with the fact that we have high 
optical depths, we do not make the distinction between the 
two sublevels. As discussed below, however, the different 
fine structure levels arc taken into accomit when choos- 
ing scattering phase functions, important for polarization 
calculations that we will present in a future paper. 

2.1.2. The scattering 

After determining the point in space where the pho- 
ton will be scattered next, we choose the thermal velocity 
components of the scattering atom. In the two directions 
perpendicular to the direction of the incident photon the 
components are drawn from a (1-D) Gaussian distribution 

with dispersion equal to ^ / . The component u„ of the 

y 

thermal velocity of the atom along the direction of the 
incident photon is drawn from the distribution 



fiVp) 



n {x — VpY + a' 



■H-^{a,x), 



(9) 



with Vp = Up{mp/2kTy/'^. To draw numbers that follow 
this distribution we use the method of Zheng & Miralda- 
Escude (2002). 

After each scattering we need to assign a new frequency 
(in the observer's frame) and direction to the photon. To 
this end we perform a Lorentz transformation of the fre- 
quency and direction of the incident photon from the ob- 
server to the atom rest frame, using the velocity of the 
atom chosen as described previously. 

Although the code ignores the level splitting with re- 
spect to the scattering cross section and the velocity dis- 
tribution, it takes into account the different phase distribu- 
tions for core versus wing scatterings, as well as for 2P1/2 
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versus 2P'^/2 scatterings. For resonant scattering, it is the 
angular momenta of the three states involved and the mul- 
tipole order of the emitted radiation that determines the 
scattering phase function. Hamilton (1940) found that 
the transition from 2Pi gives totally unpolarizcd photons 
and is characterized by an isotropic angular distribution 
function, whereas that from the 2P3/2 state corresponds 
to a maximum degree of polarization of 3/7 for a 90° scat- 
tering (also sec Chandrasekhar 1960). More specifically, 
the scattering phase function for dipole transition can be 
written as (Hamilton 1940) 



W{e) ocl + ^ cos^ B 



(10) 



with R/Q the degree of polarization for a 90° scattering 
and equal to 

(J+l)(2J + 3) 



R/Q 



(11) 



26 J2 - 15J- 1 
for the 2P3/2 ^ transition since AJ — — l,Aj = 

1, J = 3/2 according to Hamilton's conventions, and 
(2J- l)(2J + 3) 



R/Q 



(12) 



12J2 + 12J+1 
for the 2P1/2 15'i/2 transition with AJ = Aj = 0, J = 
1/2. In both equations, J is the total angular momentum 

at the excited [n = 2) state. Thus, W{6) is constant 
(isotropic) for 2Pi/2 as the excited state, whereas it equals 

H^(6') (X 1 + 3/7 cos^^ (13) 

with maximum polarization degree of 3/7 at a 90° scat- 
tering. 

On the other hand, Stenflo (1980) showed that at high 
frequency shifts (i.e., at the line wings) quantum mechani- 
cal interference between the two lines acts in such a way as 
to give a scattering behavior identical to that of a classical 
oscillator, namely pure Raylcigh scattering. Then the di- 
rection follows a dipole angular distribution with Rayleigh 
polarization 100% at 90° scattering, namely 

w{e) cx 1 + cos^ e . (14) 

Lastly, the frequency of the photon before and after scat- 
tering in the rest frame of the atom differs only by the 
recoil effect. Hence, 

(15) 



l + -^ii^(l-C0s6') 

where v, v are the frequency of the incident and scattered 
photon in the atom rest frame, respectively, the latter 
modified due to the recoil effect. This effect is negligible 
for the environments produced in the simulations. 

After determining the new direction and frequency of 
the scattered photon in the atom's rest frame we trans- 
form back to the observer's frame, and repeat the whole 
scattering procedure. 

2.2. Testing the basic scheme 

Here we present some of the tests of the RT code we 
performed against analytical solutions, as well as other nu- 
merical results that exist in the literature. In addition to 
showing the good performance of the code, these tests are 
presented here as relevant to either Ly-a emitters and/or 
the way we accelerate the code when applied in cosmolog- 
ical simulations (see §2.3). 



2.2.1. Neufeld (1990) test 

Neufeld (1990) derived an analytic solution in the limit 
of large optical depth for a source radiating resonance line 
photons in a thick, plane-parallel, isothermal semi-infinite 
slab of uniform density. The analytic emergent spectrum 
as a function of frequency shift for a midplane source is 



J{±To,x) = 



^/6 a;2 



1 



24 aro cosh[(7r4/54)0-5(|a;3 _ ifl/aro)] 

(16) 

with X = {v — vq)/ lS.V£i, Az^£) = VQ^^ksT I (mpC^) the 
thermal Doppler width, Xi the injection frequency shift 
(zero for injection at line center), and a the ratio of the 
natural to two times the thermal Doppler width. The 
quantity tq is the optical depth from midplane to one 
boundary of the slab at the line center.^ 

This analytical solution is valid in the very optically 
thick limit, with the latter being defined according to Neufeld 
(1990) as To > 10^/(v^a). This corresponds to tq > 
3.8 X 10* approximately for a temperature T=10 K as- 
sumed in the tests we present here. In deriving equation 
(16) the scattering was assumed to be isotropic. In ad- 
dition, it was assumed coherence in the rest frame of the 
atom, an assumption that makes the solution valid at the 
low density limit only, as well as approximations under 
the assiimption that wing scatterings dominate were done 
(hence the solution is valid at high optical depths). Fur- 
thermore, note that the classical slab problem is indepen- 
dent of the real size of the slab (all quantities depend on 
1/Iq with Iq the actual size of the finite dimension). Lastly, 
for this solution it is assumed that the source has unit 
strength and is isotropic, namely it emits 1 photon per 
unit time or l/47r photons per unit time and steradian. 

For center-of-linc injection frequencies the emerging spec- 
trum has maximum at x ~ ±0.88119(ay^To)"'^/^, and an 
average number of scatterings N ~ 0.909316v^to (Har- 
rington 1973, with To in these expressions defined using 
our conventions rather than Neufeld's). This scaling of 
the mean number of scatterings with optical depth in the 
case of resonant-line RT in extremely optically thick me- 
dia was first explained by Adams (1972), who understood 
that photons escape the medium after a series of excur- 
sions to the wings. Before this study it was believed that 
the number of scatterings scales with Tq , as would be pre- 
dicted by plain spatial random walk arguments (Oster- 
brock 1962). Wc briefly review the interpretation given 
by Adams (1972) with respect to the linear scaling of the 
mean number of scatterings with optical depth, since we 
refer to it extensively in the following sections. 

The mean number of scatterings is the inverse of the 
escape probability per scattering. The escape probability 
per scattering is the integral of the probability per scat- 
tering that a photon is scattered beyond certain frequency 
shift a;*. Adams (1972) identified this frequency as the fre- 

^ Neufeld's definition, used in equation (16), is such that the opti- 
cal depth at frequency shift x is given as Tx = Totj>{x). Note that 
throughout this section, with the exception of equation (16), our 
definition of tq is such that the optical depth at frequency shift x is 
given as Tx = ToH{a, x). This definition was chosen following recent 
studies (e.g. Ahn et al. 2002; Zheng & Miralda-Escude 2002) so that 
comparisons with these studies be easier. Since (t>{x) = H{a,x)/y/Tr, 
our To is smaller than Neufeld's by a factor of -v/tF. Note though that 
in the following sections we return to the Neufeld (1990) definition 

of TQ. 
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Fig. 1. — Left panel: Emergent spectra from the Monte Carlo RT {solid histograms) and as predicted analytically by Neufeld (1990) (dotted 
lines) for 3 different center— of-linc optical depths. The agreement between Monte Carlo and analytical result becomes better with increasing 
optical depth, as expected since the analytical solution is valid for very optically thick media. Right panel: The same as in left panel but in 
this case the Monte Carlo results are derived with recoil being included, whereas the analytic solution does not include recoil. The dashed 
line in the case of tq = 10'' is obtained by modifying the spectrum obtained from the Monte Carlo RT without recoil by the factor correcting 
for recoil (see text for details). 



quency where the photon, while performing an excursion 
to the wings, and before returning back to the core, trav- 
els an rms distance comparable to the size of the medium. 
Note that this is in fact an essential difference in the un- 
derstanding of resonant-line RT in extremely thick media 
compared to the spatial random walk approach. The latter 
approach assumes that during an excursion to the wings 
the photon travels an rms distance much smaller than the 
size of the medium. Thus, the first step is to determine 
X*. Using the redistribution function (i.e., the function 
that gives the probability that a photon with certain fre- 
quency shift X before scattering will have a frequency shift 
X after scattering) one can calculate both the rms fre- 
quency shift and the mean frequency shift of a photon 
which is scattered repeatedly. For a photon initially in the 
wings with a frequency shift x Osterbrock (1962) found 
that the rms shift is 1 and the mean frequency shift is 
— l/|x|. For X ^ 1, the mean shift is much smaller than 
the rms and the photon is undergoing a random walk in 



frequency with mean number of scatterings 



In real 



space, the rms distance traveled is equal to the square 
root of the mean number of scatterings times the mean 
free path. In the wings, the Voigt profile varies relatively 
slowly and the mean free path is ^ l/(/I)(.x) ^ x'^/a line 
center optical depths (we only focus on the scalings here, 
hence constants of order unity are dropped). Thus, the 
distance traveled is x/(j)(x) ^ x^ ja. Setting this rms dis- 
tance equal to tq we get .t* ^ (aro)^/^, which is in fact 
the scaling of the frequency shift where the emergent spec- 
trum takes its maximum value. Thus, going back to the 
mean number of scatterings, the escape probability per 
scattering will be ~ A{x)dx with A{x) a function to 



N. 



(t>{x)/x^dx 



be determined. According to the previous discussion x^ is 
the minimum frequency shift for which the photon during 
an excursion to the wings can travel an rms distance at 
least equal to the size of the medium. If, for simplicity, 
one assumes complete redistribution the probability that 
a photon is found after scattering with a shift between 
X and X + dx is (j){x)dx. However, this is not the prob- 
ability per scattering, since the photon will scatter ~ x"^ 
before returning to the core. Thus, A{x) is (j){x)/x^ , and 

with 4>(x) ~ a/x"^ in the wings. 

Using the above expression for one obtains Ng^ ~ tq, 
with the constant of proportionality being of the order of 
unity. 

The emerging spectra without and with recoil included 
are shown in the left and right panel of Figure 1 , respec- 
tively. A convergence test indicates that these results are 
robust if more than of order 10^ photons are used. Refer- 
ring to the left panel of the figure, the agreement between 
the results obtained with the code and the analytic so- 
lution gets better at higher optical depths. As has been 
already mentioned, the analytic solution is derived after a 
series of approximations done on the assumption of opti- 
cally thick media. For example, when deriving the analytic 
solution the Voigt function is set equal to a/irx^. Set- 
ting the Voigt function cqtial to this approximation in the 
code makes the agreement even better. The way the spec- 
trum behaves for different aro is expected qualitatively: 
the higher the optical depth or the lower the temperature 
(the higher the a), the more difficult it is for the pho- 
tons to exit the medium and the photon frequencies must 
move further away from resonance to escape. Hence, the 
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peaks of the emerging spectrum occur at higher frequency 
shifts, and the separation between the two peaks becomes 
larger. The width of the peaks gets larger with larger aro 
in agreement with the dependence of the optical depth on 
frequency (i.e., when in the less optically thick regime core 

2 

photons are relevant and the optical depth goc^s as , 
whereas in the more optically thick regime wing photons 
are more relevant, and there the optical depth scales as 
l/x2). 

In the right panel of Figure 1 we present numerical re- 
sults when recoil is included, along with the analytical 
solution (as a guide) that does not include recoil. As 
expected, including recoil shifts more photons to smaller 
(more red) frequencies. The magnitude of the effect can 
be understood as reflecting the thermalization of photons 
around frequency i^o (Wouthuyscn 1952; Field 1959). This 
process modifies the photon abundance by exp(— cc/xt), 
with xt = ksT/hAi/n- Indeed, in the right panel of Fig- 
ure 1 the dashed line for tq = 10^ is obtained by modifying 
by exp(— x/xt) the emerging spectrum obtained from the 
simulation when no recoil is included. These results arc in 
agreement with the results and interpretations by Zheng 
& Miralda-Escude (2002). 

2.2.2. Loeb & Ryhicki (1999) test 

Loeb & Rybicki (1999) address the RT problem in a 
spherically symmetric, uniform, radially expanding neu- 
tral hydrogen cloud surrounding a central point source of 
Ly-a photons. No thermal motions arc included (T = 
K). They find that the mean intensity J(f, v) as a function 
of distance from the source f and frequency shift v in the 
diffusion (high optical depth) limit is given by 

with V = vjvi,^ V = — lyphoton, the Ly-a resonance 
frequency, and i/* the frequency where the optical depth 
becomes unity. The scaled radius, r is equal to r/r*, with 
the physical distance where the frequency shift due to 
the Hubblc-likc expansion of the hydrogen cloud equals 
the frequency shift that corresponds to unit optical depth 
{= Vi,). A comparison of the results from the code with 
the analytic solution is shown in Figure 2. The analytic 
solution becomes progressively more accurate the higher 
the optical depth (or the smaller the frequency shift in the 
way this problem is parameterized, so that we are still 
in the core of the line). In addition, it deviates more 
and more from the (exact) simulation result at larger f, 
since the larger the f the more optically thin the medium 
and thus the further away we are from the assumption of 
an optically thick medium made by the analytic solution. 
Thus, the disagreement at high f is real and not an arti- 
fact caused, e.g., by small number of photons that would 
be inadequate to sample the low intensities at large f. 

2.2.3. Simple models of Ly-a emitters: Spherical clouds 
of uniform density and temperature 

Here we develop some simple models of Ly-a emitters. 
Even though there are no analytic solutions for these cases, 
one could compare our results with the published results 
of Zheng & Miralda-Escude (2002). More speficically, in 
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this section, following these authors we model spherical 
neutral hydrogen clouds. We consider two different cases 
as far as the emission is concerned. In the first case it is as- 
sumed that we have a spherical cloud with a Ly-a emitting 
point source at its center. In the second case we assume 
uniform emissivity, namely a photon is equally likely to 
be emitted from any point within the cloud. For each one 
of these two cases wc make runs assuming the cloud is 
static, contracting and expanding. In the latter two cases 
the contraction/expansion is assumed to be Hubble- like, 
namely the velocity of the neutral hydrogen atoms scales 
linearly with the radius measured from the center of the 
cloud. This velocity is set equal to 200 km/s at the edge 
of the system (and is negative/positive in the case of con- 
traction/expansion). For each case we perform two runs, 
one with column density equal to 2 x lO-'^cm"^, typical for 
Lyman limit systems, and one with column density equal 
to 2 X 10^'^cm^^. typical of Damped Ly-a systems (or line 
center optical depths equal to 8.3 x lO"' and 8.3 x 10^, re- 
spectively). In all cases the temperature is set equal to 
2 X lO^K. The initial photon frequency is assumed to be 
at the line center in the rest frame of the atom. In all 
results shown, the effect of recoil is included. Lastly, 1000 
photons were used in all runs. 

The results for the optically thin case are shown in Fig- 
ure 3 and that for the optically thick configuration arc 
shown in Figure 4. In both case the agreement with the 
results obtained by Zheng & Miralda-Escude (2002) is very 
good. These spectra can be understood qualitatively us- 
ing the way the Neufeld solution behaves depending on 
the optical thickness. In the case of an expanding cloud. 
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the photons will escape on average with a redshift be- 
cause they are doing work on the expansion of the cloud 
as they are scattered. Photons with negative frequency 
shifts (redshifted) can escape, but those with positive fre- 
quency shift (blucshiftcd) will be scattered at some point 
and they have to undergo a series of many positive shift 
scatterings to escape. Hence, the blue part of the spec- 
trum is suppressed. The situation is reversed in the case 
of a contracting cloud. It is important to keep in mind 
that the degree of suppression of one of the two peaks 
due to bulk motions depends on factors such as the op- 
tical depth and the temperature. In the case of uniform 
emissivity and expansion/contraction all spectra become 
broader because of the different velocities of the emission 
sites of the photons. In addition, when the cloud is ex- 



panding (contracting), the blue (red) part of the spectrum 
is not suppressed as much as in the central point source 
case because, at least, photons initially emitted close to 
the edge of the system have some chance of escaping even 
if they are blue (red). In the optically thicker cloud, as 
soon as the photon reaches a sufficiently large x it is not 
likely that it will be scattered by an atom with the right 
velocity to bring the photon into the line center. Rather, 
the photon will get another random shift in frequency and 
will follow an excursion in frequency while at the same 
time it diffuses spatially. This along with the fact that 
the optical depth has a power law rather than Gaussian 
dependence on x broadens the peaks compared to those of 
the optically thin case, exactly as discussed for behavior of 
the Neufeld solution. In addition, the emission peaks move 
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further away from the eenter compared to those from the 
optically thin case, since the photons have to be further 
away from the center of the line in order to escape when 
the medium is optically thick. 

2.3. Accelerating the RT 

The previous tests demonstrated that our basic Monte 
Carlo scheme works wcil for the simple test cases. When 
using it in its simplest form in high resolution cosmological 
simulations, such as the ART simulations (see §3), it takes 
unrealistic running times in order to produce results with 
sufficient numbers of photons. This is because in the case 
of resonant RT in extremely optically thick media, a sig- 
nificant amount of time is spent on the relatively insignifi- 
cant core scatterings. If we define the core through the fre- 
quency range where the Dopplcr profile dominates over the 
Lorentzian wings, then roughly speaking the core is given 

by a/-Kx1 = e~^=/\/7r, where Xc = {fc — i'o)/^^d and a 

is as defined previously. For a temperature of 10^ K, the 
core is roughly Xc = 3.5. Also, assuming complete redis- 
tribution,^ the probability per scattering for a core photon 
to exit the core is //(/ -|- erf(xc)) with / = 2a/{'KXc) and 
erf (xc) = Jq" e~*^dt. That is, roughly, the photon will 

have to scatter (7-herf(xc))/I = 1 -|-erf(xc)/I ~ 1 -M"^ ~ 
I^^ times before exiting the core. Using the above core 
definition, one finds that this is equal to ^/jTe^" / {2xc) 
and keeping only the dominant dependence on Xc, this 

2 - 

is roughly or ^ 10 scatterings. These scatterings are 
insignificant in the sense that they happen in such copious 
amounts, without being accompanied by significant spa- 
tial diffusion, since the latter occurs mostly through the 
wings. 

One way to advance photons in very high optical depths 
is to use the technique of the prejudiced first scattering 
(Cashwell & Everett 1959). With this technique one bi- 
ases the T values toward larger values than the ones that 
would be drawn from equation(l). More specifically, r is 
chosen to be uniformly distributed in [0,Tesc], with Tesc 
the optical depth for escape. Then one weights the pho- 
tons by TescC""^ to correct for the fact that r (i) is limited 
to be less than or equal to Tesc, and (ii) is assumed to be 
uniformly distributed in the < t^sc range. Using this tech- 
nique however does not improve run time requirements to 
the extent we need and clearly more drastic acceleration 
methods are needed. 

Exiting the core does not in general guarantee that the 
photons escape. In fact, the photons may return back to 
the core many times before escaping. This is not surpris- 
ing since, as we will discuss, the maximum core frequen- 
cies that can be used are much smaller than x* discussed 
previously. Especially for extremely optically thick me- 
dia (aro > 10'^), this in- and out-of-the-core procedure 
is still very expensive to follow. Hence, we accelerate our 

^ In other words, assuming that the frequency distribution after scat- 
tering is independent of the frequency before scattering and is given 
by the line profile (i.e., the source function is independent of fre- 
quency). The assumption of complete redistribution was found to 
be pretty accurate for core photons (Unno 1952; Jefleries & White 
1960). This is intuitively expected since, when in the core, the pho- 
ton frequency shift is small or comparable to the thermal velocities 
of the atoms. Thus, the latter can have a significant impa<;t on the 
frequency of the photon and in effect they redistribute it after each 
scattering according to the line profile. 



RT scheme by implementing two different methods, de- 
pending on the center-of-line optical thickness of the cell 
a photon finds itself in (tq), as well as on the thickness of 
the cell for the specific frequency shift of the incident pho- 
ton {— TQ(j}{xi) with Xi the frequency shift of the incident 
photon). In fact we parameterize the optical thickness of 
a cell not only via tq, the line-of-center optical depth from 
the center of the cell to one of its edges, but rather via the 
product of a and tq, motivated by the Neufeld solution. 
This parameterization turns out to be very good for media 
less optically thick than those the Neufeld solution applies 
to. We discuss these two acceleration methods, as well as 
some additional acceleration techniques in the following 
subsections. 

2.3.1. Extremely optically thick cells: Controlled Monte 
Carlo motivated by the Neufeld solution 

This acceleration scheme is based on controlled Monte 
Carlo simulations of resonance RT in cells (cubes) with 
several physical conditions, representative of the extremely 
optically thick cells in the simulations. The idea is to ob- 
tain trends and best-fit functional forms for the spectra 
emerging from thick cells. These spectra can then be used 
when running the code so that instead of following the 
scattering of the photons in detail, we can draw the fre- 
quency of the photon emerging from a thick cell using the 
pre-calculated spectrum appropriate for the physical con- 
ditions in this cell. In principle, controlled Monte Carlo 
simulations can be used for any range of optical thick- 
nesses. We use it only in the extremely optically thick 
cells where (aTo)e// > 2 x 10^ and {TQ(f){xi))ef j 1 with 
Xi the frequency shift of the incident photon.^ We do that 
because we are motivating this method by the Neufeld so- 
lution which is applicable only at the diffusion limit. The 
inherent cell structure of the AMR simulation outputs or 
the cell structure that can be generated for other cosmolog- 
ical codes, along with the resolution imposed isothermality 
and uniformity of each cell, are conducive to some kind of 
modification of the Neufeld solution. In some sense, with 
the advent of cosmological simulations, the contemporary 
analogue of the extensively studied classical slab problem 
is the completely unexplored problem of resonance RT in 
a (;ube. This motivated a detailed study of the resonance 
RT problem in cubes where the reader is referred to for 
more details and results (Tasitsiomi 2006). Here we only 
summarize briefly some key results relevant to the current 
study. 

As discussed in §2.2.1, the Neufeld solution was obtained 
under some assumptions. To fit the controlled Monte 
Carlo spectra with a Neufeld type spectrum we have to in- 
vestigate how sensitively the solution depends on these as- 
sumptions, as well as whether these assumptions arc valid 
in cosmological simulations. This is done in the following 
paragraphs. 

^ We have used the index eff because, as is discussed later in this 
section having in mind an implementation of the RT code for AMR 
simulations, to decide whether this method is applicable or not we 
create a mesh on top of the simulation mesh. In this new mesh, the 
photon is always at the center of a cell. Then it is the 'effective' 
physical conditions in this new cell that are relevant when decid- 
ing if the acceleration method at hand is applicable or not. In the 
case of simulations without a cell structure, the index eff becomes 
redundant, since there is no initial mesh to begin with. 
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Choosing the exiting frequency 

The exiting frequency of a photon entering an extremely 
optic;ally thick cell is drawn by an emerging frequency dis- 
tribution similar to the Neufeld solution (equation 16). 
However, the Neufeld solution is derived for a semi-infinite 
slab, whereas the simulation cells are finite cubes. Further- 
more, the solution assumes isotropic scattering, no recoil, 
which anyway is negligible in the simulations, and does 
not include velocities such as those associated with pecu- 
liar motions or the Hubble expansion. Lastly, it assumes 
that the source of the radiation lies within the slab,^ and 
is valid for optically thick frequencies {To(j){xi) » 1). 

Starting from the point on bulk velocities, we use the 
Neufeld solution - applicable for an observer moving with 
the bulk flow of the fluid - by taking into account the 
way the specific intensity transforms between two incr- 
tial observers moving at a certain speed with respect to 
each other (i.e., Iv/v'^ is invariant, where Iv is the number 
of photons rather than the energy intensity. In the latter 
case, the quantity that would be invariant would be J^/^'') . 
The second point wc address has to to do with the slab 
versus cube difference between the analytical solution and 
the simulations. As discussed in §2.2.1, Neufeld's solution 
depends on one parameter, aTQ. Qualitatively, one expects 
that the spectrum emerging from a cube rather than a slab 
be well described by the same solution but for an effective 
aro smaller than the actual aro of the cell. The reason for 
this is that when for example observing the emergent flux 
from the z-dircction in a cube, wc lose all photons that 
in the case of the slab would wander, scatter many times 
along the infinite dimensions and finally find their way out 
from the z-plane. In the case of the cube these photons 
would not be counted simply because they have exited the 
cube from planes other than the z-plane. This would be 
equivalent to solving the problem that Neufeld solved but 
this time including losses of photons (or, more appropri- 
ately, by generalizing the 2-dimensional diffusion equation 
derived by Neufeld into a four dimensional one - instead 
of r, v now the intensity will be a function of Tx,Ty, and 
p). Numerical experimentation of RT in cubes and slabs of 
the same physical conditions, verified that the above guess 
is correct. In fact, the cube spectrum is well described by 
the Neufeld solution for a slab if 2/3 of the aro of the cube 
are used as input parameter to the slab analytic solution. 
This is shown in Figure 5 (also see Tasitsiomi 2006). 

Furthermore, the Neufeld solution assumes that the source 
of radiation lies within the slab (or cube in our case). In 
fact, the version of the solution we have been discussing so 
far (equation 16) assumes that the source is at the center 
of the slab. However, in the case of mesh-based codes, as 
photons cross from one cell to the other, in general the 
source is not at the center of the cell. For codes without 
an inherent cell structure the obvious solution to this is 
to create a cell and have the photon at each instant at 
the center of the cell. As is discussed in what follows, this 
turns out to be the most efficient solution in the case of 

^ More specifically, equation 16 assumes that the source is a plane 
source in the middle of the slab. Due to symmetry arguments, a 
plane source located at the middle of a slab is equivalent with respect 
to the spectrum of the emergent radiation to a central point source. 
Neufeld (1990) provides a more general expression for different source 
positions. 
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mesh-based codes as well. 

Neufeld provides a more general expression for various 
source positions within the slab, as well as for the trans- 
mission and reflection coefficients assuming an external 
source. Using either option for mesh-based codes, trying 
to take advantage of the already existing mesh structure, 
creates complications: in the case of a non-central but in- 
ternal source, the equivalence of a point or infinite plane 
source - necessary for all the above discussion to be valid - 
breaks down if the source is not located at the center of the 
slab. And using the reflection/transmission probabilities 
makes the algorithm more complicated. But most impor- 
tantly, there is an intrinsic limitation in the simulations 
due to finite resolution: it is not clear how meaningful it 
is to be discussing differences in position less than the cell 
size (i.e., if one can really tell the edge from the center of 
the cube). Instead, at every point the photon is found we 
create a new mesh on top of the simulation mesh. The pho- 
ton is always found at the center of a cell whose physical 
parameters are calculated using the cloud-in-cell weighting 
scheme. Each time the size of the cell is set to the simula- 
tion cell size the photon is in. Note that it is the physical 
parameters of this effective cell that determine the way the 
code proceeds (i.e., if the effective cell aro is larger than 
2 X lO"^ and To0(a;j) :» 1 then the controlled Monte Carlo 
results are used. If one of these two conditions (or both) 
is not satisfied in the effective cell then the code returns 
to the original cell. Depending on the original cell phys- 
ical conditions and the photon frequency either the exact 
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Monte Carlo or the method described in §2.3.2 is used). 

In the Neufeld solution the condition To(f>{xi) ^ 1 allows 
him to truncate a series appearing in the solution process 
by keeping up to first order terms in l/{To(f>(xi)). Thus, 
the solution is valid only for optically thick injection fre- 
quencies. We find that the higher order corrections are 
pretty small. However, for a certain tolerance, one must 
decide how thick is thick enough for the Neufeld spectrum 
to be applicable. We take that the spectrum from a slab is 
satisfactorily predicted by the analytical solution for fre- 
quency shifts for which ro^(xj) > 10. 

As has been shown in §2.2.1 the recoil effect can be eas- 
ily accounted for multiplying the Neufeld solution by the 
appropriate factor. In any event, the recoil effect for our 
conditions is negligible and hence is dropped in the simula- 
tion calculations. To see this, the recoil effect corresponds 
to a frequency shift that would be caused by a velocity 
~ hv/m,pC = 3 m/s. This velocity is negligible compared 
to the thermal velocities expected in cosmological simula- 
tions, and given the peculiar and Hubble flow velocities, 
the small non-coherence in the atom's rest frame intro- 
duced by recoil will be totally unobservable. Hence, the 
Neufeld approximation is good in that respect as well. 

Choosing the exiting direction and point 

Referring to /x, the cosine of the angle with which the pho- 
ton is exiting a cell, measured with respect to the normal 
to the exiting surface, we draw its value from the following 
cumulative probability distribution function (cpdf ) (Tasit- 
siomi 2006) 

P(</.) = ^(3 + 4m). (18) 

This cpdf is found to be an excellent description of the 
directionality of the emergent spectrum and clearly devi- 
ates from isotropy. In fact, it verifies the findings of other 
studies that in optically thick media photons tend to exit 
in directions perpendicular to the exiting surface (see, e.g., 
Chandrasekhar 1960; Phillips & Mcszaros 1986; Ahn et al. 
2002). In the case of RT in accretion disks this has been 
identified as an expected limb darkening (or 'beaming'; 
Phillips & Meszaros 1986) of the disk (i.e., the disk is very 
bright when observed face on and less bright when ob- 
served edge on). In cases of very optically thick media, the 
emerging radiation directionality approaches the Thomson 
scattered radiation emergent from a Thomson-thick elec- 
tron medium. This Thomson limit obtained initially by 
Chandrasekhar (1960), was confirmed later numerically by 
Phillips & Meszaros (1986). 

It has been implied by some authors (Ahn et al. 2002) 
that the fact that in optically thick media RT occurs mostly 
via wing photons with the latter being described by a 
dipole phase function (see §2.1.2), and the fact that Thom- 
son scattering is also described by a Rayleigh (dipole) 
scattering phase function, explains why the resulting /x 
probability distributions are similar. However, we find the 
same cpdf when the scattering is taken to follow either an 
isotropic or a dipole distribution. For such optical thick- 
nesses the details of the exact phase function do not mat- 
ter, at least not with respect to the exiting angle cpdf. All 
the phase functions involved in Ly-a scattering are only 
mildly anisotropic and they simply enhance a little bit the 



coherence of the scattering at the observer's frame com- 
pared to the isotropic scattering case. So the fact that 
the exiting angle cpdf in extremely optically thick slabs 
(cubes) does not depend crucially on the assumptions on 
the phase functions docs not come as a surprise. The un- 
derlying physics is simply that in extremely thick media 
most of the photons escape along the normal to the slab 
where the opacity is smaller. The azimuthal angle (p with 
which the photon exits a cell is distributed fairly uniformly 
in [0, 27r] (for more details see Tasitsiorni 2006). 

Referring to the distribution of exit points, one can ar- 
gue that trying to specify the exact coordinates of the exit 
point of a photon from a simulation cell is, in some sense, 
superfluous since there is always the resolution limitation. 
Thus, we assume that the exiting points are distributed 
uniformly. The deviations of the exiting points from uni- 
formity are relatively small (Tasitsiomi 2006). Similarly, 
resolution limitations make us focus on total distribution 
functions of photon properties - where total here means 
distributions averaged over an entire cube side - without 
regards to a possible dependence of these distribution func- 
tions on the photon exit point. 

Lastly, we have checked whether the emergent photon 
parameters can be drawn independently. We found no 
significant correlations among them (e.g., we checked for 
correlations between emergent frequency shift and (pre- 
ferred) range of exiting directions). Thus, drawing them 
independently is correct. 

2.3.2. Moderately optically thick cells: Skipping the core 

scatterings 

This acceleration scheme is used if the cell the photon 
is in has 1 < aro < 2 x 10^. It is also used in the case 
of cosmological simulation codes with a pre-existing mesh 
when the cell the photon is in has aro > 2 x 10*^, but 
the effective cell (see §2.3.1) has 1 < aro < 2 x 10"^, and 
thus the previous acceleration scheme (discussed in §2.3.1) 
is not applicable. The scheme is based on the idea that 
if a photon is within a certain core (to be determined), 
we can skip all the core scatterings and go directly to the 
scattering with a rapidly moving atom that can bring the 
photon out of the core (for some first implementations of 
this idea see Avery & House 1968; Ahn et al. 2002). As 
soon as this happens, the initial detailed transfer resmnes 
until either the photon escapes or re-enters the core. The 
scheme's validity relies upon the correct choice of the core 
value, so that while in the core the photon does not diffuse 
significantly in space, whereas significant diffusion occurs 
when the photon exits the core. 

To achieve the scattering that brings the photon out- 
side the core we choose thermal velocities (in units of 
\j2kT/m) from the distribution (Avery & House 1968; 
Ahn et al. 2002) 

p{v) = ^e-^' (19) 

and in the range [vmimi'max]- The lower limit Vmin is 
the minimum velocity necessary for the photon to just 
make it to the core Xc- The upper limit is formally in- 
finite, but for any practical realization it can be set to a 
large enough number (e.g., \/ x1 + 10). For a scattering to 
bring the photon to just Xc from the center, independent 
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of the directions of incident and outgoing photon, and un- 
der the assumptions of colierence in the rest frame of the 
atom, isotropic scattering phase function, and zero radia- 
tion damping, it can be shown that Vmin = TOaa::(|a;|, \xc\) 
(Hummer 1962), with x the initial frequency shift (as usual 
in units of the thermal Doppler width). In our case it is 
always Vmin = \xc\ since the photon is inside the core. 
We checked and verified that the assumptions under which 
Vmin is derived are good for cosmological simulations. This 
is not surprising since, e.g., the assmnption of an isotropic 
phase function is not very crucial. As discussed already, 
none of the relevant phase functions is strongly anisotropic. 
Those that are anisotropic simply tend to favor slightly 
smaller frequency shifts (since they favor post-scattering 
directions close to pre-scattering directions) and hence in- 
crease a little bit the coherence in the observer's frame 
from scattering to scattering. At the limit of many scat- 
terings (and while still at the optically thick regime) this 
is not a significant effect (for the tiny differences in the fre- 
quency redistribution function with isotropic versus dipole 
phase function see Figure I of Hummer 1962). Or, the as- 
sumption of coherence in the rest frame of the atom is 
also expected to be a pretty good assumption for the me- 
dia in the simulations from the point of view of the recoil 
effect, as we discuss in §2.3.1, and from the point of view 
of collisions as we discuss in §3.5.1. 

To motivate the core values we can use (i.e., the maxi- 
mum frequency shifts for which we can ignore the repeated 
scatterings without biasing the results) we must take into 
account the different physics of resonant RT in the two 
different regimes, 1 < utq < 2 x 10^ and utq > 2 x 10^. In 
the first regime photons escape on a single longest flight 
(Adams 1972) in accordance with the understanding of 
resonant RT in moderately thick media developed by Os- 
terbrock (1962). In this thickness regime the important 
frequency is the frequency where the optical depth be- 
comes unity. Photons within this frequency shift barely 
diffuse in space, whereas as soon as they exit this fre- 
quency shift they escape while taking their longest spa- 
tial step {flight). In the second, extremely optically thick 
regime (aro > 2 x 10"^) as Adams (1972) suggested, pho- 
tons escape during a single longest excursion rather than 
flight. In this case the important frequency is the fre- 
quency with the following property: if a photon is given 
this frequency and is left to slowly return to the center of 
the line (by performing a double random walk, in space 
and frequency) , the overall rms distance that it will travel 
in real space while returning to the line center equals the 
size of the medium (i.e., the important frequency shift in 
this case is the shift discussed in §2.2.1). This physics 
motivates our cores, i.e., for moderately thick media the 
core must be safely optically thick, whereas for extremely 
optically thick media the core must be safely smaller than 
X*. Then using numerical experimentation we find the ex- 
act maximum possible core values that can be used in each 
case. 

A comparison of the exact Monte Carlo and the core 
acceleration scheme applied to moderately thick media is 
shown in Figure 6. Note that these spectra are one cell 
runs, and are not the final results of the RT around the Ly- 
a emitter (which are discussed in a later section). In the 
top panel, we present the exact emergent spectrum from 
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Fig. 6. Top pane/.- Comparison of the exact Monte Carlo results 
('ex.', solid line) and the results obtained using the eore acceleration 
method ('ap.', dotted line) for the minimum cell qt = 1 for which 
this acceleration method is used in the Ly-a RT code. Also shown is 
a larger core frequency, Xc = 0.2, which shows the way the emergent 
spectrum is biased if one uses a higher core frequency. Bottom panel: 
Same as in top panel but for more optically thick cells, aro = lO'^. 
In this case a core frequency Xc = 0.8 can be used. In addition, 
we show exact results for a different pair of temperature and optical 
depth {dashed line), that however correspond to aro = 10^. Clearly, 
aTQ is a good way to parameterize the problem at these moderate 
optical thicknesses. 



a cube with aro = 1, as well as the spectrum obtained if 
a core Xc = 0.02 is used. Despite it being a pretty small 
core, it improves the speed of the algorithm by orders of 
magnitude.* Also shown is what the bias would be if one 
used a higher core frequency {xc = 0.2): photons would 
be artificially shifted at higher (absolute) frequency shifts. 

^ The exact improvement factor depends on optical thickness, and 
is higher for thinner cells. Furthermore, the improvement factor is 
different for the same aro but different temperatures and optical 
depths. More specifically, it is higher for lower optical depths and 
temperatures. 
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To find the maximum core that can be used without this 
biasing, we made runs with successively higher cores. We 
use as cores: 0.02 for 1 < aro < 10, 0.1 for 10 < aro < 10^ 
and 0.8 for 10^ < aro < 2 x 10^. One can easily verify 
that for a wide temperature range these cores are safely 
within the optically thick regime. 

The comparison between the exact Monte Carlo and 
the accelerated scheme for optically thicker cells (but still 
at the moderately thick regime) is shown at the bottom 
panel of Figure 6. We have seen via the Neufeld solution 
that characterizing a slab - or a cube in our case - us- 
ing aro is very good in the case of very optically thick 
media (aro ^ 10^). In the bottom panel of Figure 6 we 
present two different sets of temperature and tq, which 
nevertheless correspond to the same aro (and smaller than 
that for which the Neufeld solution is applicable). Clearly, 
aro parameterizes nicely enough these emergent spectra 
as well. This fact justifies our classification of simulation 
cells with respect to their aro value. Note that the fact 
that the emergent spectrum for these physical conditions 
seems to depend on aro is not trivial, and was checked 
only for ranges of temperature and optical depth that are 
anticipated to be relevant to cosmological simulation en- 
vironments. A simple way to see why this may not be a 
general statement comes from the physics of RT in mod- 
erately thick media. As discussed in such media photons 
escape roughly when they reach the frequency where the 
optical depth is unity. If, for example, the frc;quency shift 
X where the optical depth becomes unity is within the 
Doppler core (as anticipated) then this frequency shift is 
defined through roe^^ = 1 and clearly depends only on 
ro and not on temperature. This is in contrast to ex- 
tremely optically thick media where the frequency shift 
relevant for escape through the single longest excursion is 
.T* ^ {aToY^^ (see §2.2.1), namely it depends on aro. 

In the case of extremely thick media we find roughly 
the following maximum possible cores: 3 for 2 x 10^ < 
aro < 10*, 5 for 10* < aro < 10^ 7 for 10^ < aro < 10^ 
17 for 10^ < aro < 10^, 30 for 10^ < aro < 10^ and 
80 for aro ^ 10^- As an example, in Figure 7 we show 
the Neufeld prediction for the emergent spectrum from a 
slab with aro = 10'' and the results of our acceleration 
scheme using a core Xc = 30. This is a quite large core 
frequency, and still the acceleration scheme gives a very 
accurate emergent spectrum. The core values we find scale 
with x■^, roughly as Xc — 0.15a;*. 

The aro-dependent core frequencies that we motivate 
here based on the different physics for different aro regimes 
is a quite new approach. Previous studies (e.g., Hansen & 
Oh 2005) define the core frequency as the frequency where 
the wings start dominating over the Doppler core. Clearly, 
to achieve the best efficiency of the acceleration scheme, 
which is highly desirable in our applications due to the very 
complex environments, we have to use a depth-dependent 
core definition. Other authors who considered variation 
of the core frequency with temperature and optical depth 
(Ahn et al. 2002) find a bit different values than ours, at 
least for the low aro range that they worked with: they 
find that a core frequency of about ^/^^ can be used for 
aro > 10^, with slightly higher values permitted for even 
larger ro- However, we find that this value is a bit large for 
aro — 10^, and that significantly higher core values can be 
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Fig. 7. — Comparison of the analytic solution obtained by Neufeld 
(1990, 'ex.', solid line) and the results obtained using the approxi- 
mate core acceleration method ('ap.', dotted line) for ar = 10'^ and 
Xc = 30. The point of the figure is that at extremely high optical 
depths the core values one can use can be pretty high. 



used for higher ro . The reasons for our disagreement with 
Ahn et al. (2002) are not clear. 

The discussion with respect to the validity of this ac- 
celeration scheme has been limited so far to the emergent 
spectrum of radiation. One would expect to see that in- 
deed the assumption that the photons do not move signif- 
icantly in space during the multiple core scatterings that 
are skipped is true. And, that all other quantities, such as 
exit point and exit angle distributions remain the same, in 
addition to the emergent spectrum. The latter has been 
tested and found true. Furthermore, note that the an- 
gle information is relevant mostly when the photon is at 
the optically thin regime, where anyway we use the exact 
transfer scheme. With respect to the exit points, or dis- 
tances that the photons move while in the core, since these 
are not larger than one cell size, limitations due to the fi- 
nite simulation resolution render these concerns moot. To 
get an idea, following an argument similar to that pre- 
sented in §2.2.1 leading to ~ (aro)^/^, and using the 
scaling Xc — O.lSx* one finds that by ignoring the scat- 
terings within the core for extremely optically thick cells 
roughly one ignores a spatial diffusion of the photons of 
order 10^"^ of the size of a simulation cell. 

In summary, each time a photon enters a simulation cell, 
there are the following three possibilities: 

1. If the cell has aro < 1) the exact Monte Carlo RT 
is used. 

2. If the cell has 1 < aro < 2 x lO'^ and the photon 
frequency shift is |a;| < Xc, then we skip the core 
scatterings. If the photon frequency is outside the 
core we use again the exact Monte Carlo RT. 
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3. If the cell has aro > 2 x 10"^ then if there is no pre- 
existing mesh structure of the cosmological simula- 
tion then if the frequency of the photon is such that 
ro0(x) » 1 then the controlled Monte Carlo results 
arc used. If there is a pre-existing mesh structure, 
then if aro > 2 x 10"^ then the physical conditions 
of the effective cell are calculated. If for the effec- 
tive cell it is (i) aro > 2 x 10'^ and the frequency 
of the photon is such that (ii) To0(a;) » 1, then we 
use the controlled Monte Carlo motivated by the 
Neufeld solution. If either (i) or (ii) is not true, 
then the first acceleration scheme is tried for the 
original rather than the effective cell. And if it is 
not applicable, then the exact Monte Carlo scheme 
is used. 

2.3.3. Calculating images and spectra 

To construct images of the Ly-a emitters for various di- 
rections of observation the code calculates the contribution 
to the image along a certain direction at each scattering 
(see, e.g. Yuscf-Zadeh et al. 1984; Zheng & Miralda-Escudc 
2002). This contribution is e~^'^'= P^tj), ji) where t^sc is the 
optical depth for escape from the current scattering po- 
sition along the direction of observation to the observer, 
/i is the cosine of the angle between the direction of the 
incident photon and the direction of observation, (p is the 
azimuthal angle, and P{4>i At) is the normalized probability 
distribution for the photon direction (in fact P is indepen- 
dent of in our case). 

This way of calculating images and spectra has the ad- 
vantage of giving fairly good statistics for relatively small 
numbers of photons. Thus, by lowering the number of 
photons needed for the results to converge, it can poten- 
tially speed up the calculations. It also converges rapidly 
for the fainter parts of the source, hence it is very useful 
for sources with high emissivity contrast. One disadvan- 
tage is that due to computing resources limitations it lim- 
its the calculations to only a small number of pre-chosen 
directions of observations. In addition, for complicated ge- 
ometries such as those produced in simulations one must 
verify that running more photons is more expensive than 
calculating t^sc used in this method. We find that indeed 
this is the case for the ART environments where the RT 
code is applied in this study. 

2.3.4. Parallelization 

To reach high performance we implement the parallel 
execution of the code. Our Monte Carlo scheme is partic- 
ularly easy to parallelize, since each ray is independent of 
others. The parallelization is done using the Message Pass- 
ing Interface (MPI) library of routines. As every photon 
ray is independent, communication requirements among 
the different processes are minimal, and in essence MPI 
distributes copies of the code which are run autonomously 
in the different nodes used. However, each processor is as- 
signed and runs photons from different emission regions. 
To get an idea about the performance of the code (using 
the above acceleration schemes), 10^ photons ^ transfer to 
10 physical kpc from the center of the ART Ly-a emitter 

^ This number of photons is well above the minimum necessary for 
the results to converge as will be discussed in a later section 



we apply the code to in about 4 hours on 8 Intel Xeon 3.2 
GHz processors on the Tungsten NCSA cluster . 

2.4. Final images and spectra of simulated Ly-a emitters 

The detailed Ly-a RT is carried out up to a certain 
distance from the center of the source and then the Ly-a 
GP absorption is added. This distance where the detailed 
RT stops is determined through a convergence test. The 
existence of such a scale is guaranteed given that the fur- 
ther away a photon moves from the center of the object, 
the most improbable it becomes for it to scatter back in 
the direction of observation. Furthermore, the size of this 
convergence radius can also be motivated observationally, 
from the extent of Ly-a halos that have been observed. 

The surface brightness of each pixel of the constructed 
image is 
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a 



X e 



(20) 



where the sum is over the fluxes of all photons (i), and all 
their scatterings (j) with scattering positions that project 
onto the pixel; ^pix is the angle subtended by the pixel to 
the observer, and the factor e""^"^^ accounts for the dimin- 
ishing of the brightness due to the hydrogen intervening 
between the radius where the detailed RT stops and the 
observer. To find the flux Fi_j carried by each photon at 
each interaction, we first calculate the total luminosity, 
Ltot, of the emitter through the sum of the luminosities of 
the individual source cells. For N photons (or more ac- 
curately wavepackets) used in the Monte Carlo, then each 
photon carries a luminosity Fi_j (independent of photon 
and scattering numbers i and j, respectively, in our case) 
equal to 

= ¥4 

where dL is the luminosity distance calculated for the 
adopted cosmology. Note that there is no l/47r factor. 
This factor comes from P{4>, ji) - in equation (20) - which 
is normalized to unity. 

The GP absorption optical depth is calculated as de- 
scribed in Hui et al. (1997). It is calculated for each pixel 
separately, and the number of different lines of sight that 
have to be used per pixel is determined by checking con- 
vergence of the final result. For high enough image spatial 
resolution (similar to the one used in this study) one line 
of sight per pixel is enough, since the simulations them- 
selves have finite spatial resolution. The characteristics of 
the line emerging after the detailed RT (i.e., its width) and 
before adding the GP absorption determine how far away 
in distance one must go when calculating tqp, since one 
needs to go up to the point where the shortest line wave- 
length is redshifted at least to the Ly-a resonance because 
of Hubble expansion. Often, this physical distance is larger 
than the physical size of the cosmological simulation box. 
In this case, we take advantage of the periodic boundary 
conditions and use replicas of the same box making sure 
we do not go through the same structures. This turns out 
to be easily done as long as one does not have to use the 
box too many times (more than ~ 5). Furthermore, we 
consider two distinct scenarios, one where the effect of the 
red damping wing is taken into account and one where the 
red damping wing is suppressed as would be the case if for 
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example the Ly-a emitter was in the vicinity of a bright 
quasar. 

Lastly, spectra are obtained by collapsing the 3-D image 
array (2 spatial dimensions+wavelength) along the spatial 
dimensions. 

3. APPLICATION TO COSMOLOGICAL SIMULATIONS 

3.1. The simulations 

Here we present some basic information regarding the 
cosmological simulations we use in what follows in order 
to apply the Ly-a RT code in a cosmological setting. 

The RT is carried out using outputs of the ART code for 
the concordance flat ACDM model: CIq = 1 — fl\ = 0.3, 
h = 0.7, where VIq and f^A are the present-day matter and 
vacuum densities, and h is the dimcnsionlcss Hubble con- 
stant defined as i?o = lOO/i km s^^ Mpc^^. For the power 
spectrum normalization the value ag — 0.9 is used. This 
model is consistent with recent observational constraints 
(e.g., Spcrgcl et al. 2003). The initial conditions of these 
simulations are the same as those in Kravtsov (2003) and 
Kravtsov & Gnedin (2005), loading to the formation of a 
Milky Way sized galaxy at 2; = 0. However, these simu- 
lations are different in that, in addition to dark matter, 
gas dynamics, star formation and feedback, cooling, etc., 
they also include non-equilibrium ionization and thermal 
balance of H, He, H2 and primordial chemistry, full RT 
of ionizing radiation and optically thin line RT of Lyman- 
Werner radiation. The continuum RT is modeled accord- 
ing to the Optically Thin Variable Eddington Tensor ap- 
proximation described in Gnedin & Abel (2001), whereas 
cooling uses the abundances of species from the reaction 
network, as well as corrections for cooling enhancement 
due to metals. 

The code reaches high force resolution by refining all 
high-density regions with an automated refinement algo- 
rithm. The criterion for refinement is the mass of dark 
matter particles and gas per cell. Overall there are 9 re- 
finement levels. The physical size of a cell of refinement 
level I is 26.161 x 2^~' pes at z ~ 8 (the rcdshift wo fo- 
cus on in this study). The dark matter particle mass at 
the highest resolution region is 9.18 x lO^/i~^M0, and the 
box size for which results are presented in this paper is 
6ft.-iMpc. 

For each simulation cell we have available information 
such as the temperature, the peculiar velocity, the neutral 
hydrogen density, the ionized hydrogen density, the metal- 
licity, etc. With this information and using the mesh of the 
ART code itself we follow how Ly-a photons are being ini- 
tially emitted and subsequently getting scattered. As an 
example of an application of the Ly-a RT code developed 
for the ART code we focus on the most massive emitter 
at z ~ 8. This emitter is found within a highly ionized, 
butterfiy-shaped bubble. Outside this bubble the Uni- 
verse is highly neutral, whereas some dense neutral cores 
associated with the forming galaxy exist within the bub- 
ble. Results for more emitters, different redshifts, multiple 
directions of observation, larger simulation boxes, etc., will 
be presented in future papers. 

3.2. Intrinsic Ly-a emission 

There are a number of different mechanisms that can 
produce Ly-a emission from high-redshift objects. Here 



we classify them into recombination and collisional emis- 
sion mechanisms. By recombination emission mechanisms 
we refer to Ly-a photons that are the final result of the cas- 
cading of recombination photons produced in ionized gas. 
The gas may be ionized by the UV radiation of hot, young, 
massive stars, from an AGN hosted by the galaxy, or by 
the intergalactic UV background. By collisional emission 
mechanisms we refer to photons that are produced by the 
radiative decay of excited bound (neutral) hydrogen states, 
with collisions being the mechanism by which these excited 
states are being populated. This mechanism takes place 
when gas within a dark matter halo is cooling and col- 
lapsing to form a galaxy and radiates some of the gravita- 
tional collapse energy by coUisionally excited Ly-a emis- 
sion, when gas is shock heated by galactic winds or by 
jets in radio galaxies, and in supernova remnant cooling 
shells. We underscore the fact that the states are bound 
states, because in principle collisions can also cause ion- 
ization in which case we would have production of Ly-a 
photons under a recombination mechanism, according to 
our definition conventions. With the exception of AGN 
and jets, which are not included in ART simulations, as 
well as the fluorescence emission due to the intergalactic 
UV background which would be relevant at lower redshifts 
than we focus on in this study, we will try to briefly assess 
the importance of these separate Ly-a emission sources. 
This is interesting in particular because, in addition to 
the different dependenc;e on the physical parameters (i.e., 
different temperature dependence and dependence on ion- 
ized versus neutral hydrogen), these mechanisms may also 
have a different spatial distribution. For example, shock 
heated gas from gravitational collapse may be a spatially 
more extended Ly-a source than the gas photoionized by 
UV radiation of young stars at the relatively compact star 
forming regions. The dominant source of Ly-a emission 
may be what distinguishes most Ly-a emitters from the 
more extended sources referred to in literature as Ly-a 
blobs (Steidel et al. 2000; Haiman et al. 2000; Fardal et al. 
2001; Bower et al. 2004). 

Before discussing the different Ly-a emission mecha- 
nisms, we should first mention that, due to practical lim- 
itations (i.e., we can only use a relatively limited number 
of photons), we use as source cells only the cells that con- 
tribute significantly to the total luminosity of the object. 
Hence, we set a threshold on the cell luminosity and use 
as source cells only the cells whose luminosity exceeds this 
threshold. Then by performing a convergence test, namely 
by doing runs assuming different luminosity thresholds up 
to the point where including lower luminosity source cells 
does not change the results (within some pre-specificd tol- 
erance), we determine the minimum luminosity a simula- 
tion cell must emit to be one of the cells where photons will 
originate from. It is meaningful to consider a similar con- 
vergence check with respect to the Ly-a RT results, and 
this will be discussed in a later section. The convergence 
test reveals that the luminosity of the object is dominated 
by a few very luminous cells. To get an idea, the luminosi- 
ties of cells within the virial extent roughly range from 
10'*"'^ to several times 10"'''* photons/s. The total luminos- 
ity of the object is the sum of the luminosities of the cells 
considered. Even though most of the volume, say, within 
the virial radius is in low to moderate luminosity cells. 
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the sum of the luminosities of these cells is not significant 
enough compared to the less numerous high luminosity 
cells. For the object at hand the convergence test suggests 
that one can use as source cells only cells with luminosities 
above ~ 5 X 10'''" photons s~^. This value determines the 
relative importance of the different Ly-a emission mecha- 
nisms discussed in what follows. With the aforementioned 
luminosity threshold, the total luminosity of the emitter at 
hand is roughly equal to 4.8 x 10^^ ergs/s. We sample the 
emission rcigion (i.e., the cells with luminosity al)ovc the 
luminosity threshold discussed) by emitting equal weight 
wave packets, but in numbers that reflect the relative lu- 
minosities of the cells. 

Note that this discussion on the various mechanisms, 
emission rates, etc., should somehow be affected by the 
limited simulation resolution, a factor that will be studied 
in detail in the future. Furthermore, the approach adopted 
in this section is an 'ordcr-of-magnitude' one. We defer a 
more thorough and statistical analysis of the Ly-a emis- 
sion sources in high redshift galaxies to a future study, 
where all factors will be taken into account. For exam- 
ple, the discussion about the importance of the various 
emission mechanisms must be extended to the after RT 
results and after including dust. This is because it could, 
for example, be the case that recombination Ly-a photons, 
despite being more numerous as discussed below, may be 
more likely to be absorbed than collisional Ly-a photons, 
if one assumes that there is more dust in star forming re- 
gions - where recombination photons are generated - than 
in regions where collisional Ly-a photons originate from. 

3.2.1. Ly-a photons from recombinations 

The recombination rate of a cell is 

r = neUpaBV (22) 

with UejUp the number density of electrons and protons, 
respectively, and V the volume. In principle species other 
than hydrogen may contribute to Ue- Thus, rig in general 
is not equal to rip. In what follows, we take into account 
electrons contributed by the ionization of He. Other BBN 
predicted species such as Li, Be and B (with, anyway, tiny 
abundances), and elements produced through stellar pro- 
cessing such as C, N and O are not taken into account. 

Recombination photons are converted with certain effi- 
ciency into Ly-a photons. In particular, for a broad range 
of temperatures centered on T = 10'* K, roughly 38% of 
recombinations go directly to the groimd state. A fraction 
^ 1/3 (32%) of the recombinations that do not go to the 
ground state go to 25* rather than 2P and then go to the 
ground state via two continuum photon decay (cf. Table 
9.1 of Spitzer 1978). Hence, only a fraction ^ 40% of the 
recombinations yield a Ly-a photon. The temperatures 
of simulation cells within the virial extent of the emitter 
are in the 10^-^ — 10^-^ K range, with most cells in the 
10* — 10^ K range. Due to the weak temperature depen- 
dence of the various recombination coefficients the above 
conversion efficiencies are roughly applicable throughout 
this temperature range. Furthermore, if the gas is optically 
thick, then photons that originate from recombinations to 
the ground state will be immediately absorbed by another 
neutral hydrogen atom and eventually they, as well, will 
produce Ly-a photons. Assuming for now that this is the 
case (as will be discussed later in this section), as well 



as that the medium is thick in Lyman-series photons, so 
that all higher Lyman-series photons are re-captured and 
eventually yield Ly-a photons, we adopt case B recombi- 
nation. For the recombination coefficient we use the fit 
obtained by Hui & Gncdin (1997), accurate to 0.7% for 
temperatures from 1 to lO'^K 

as = 2.753 x lO-^cm^s"! ^' ',,,22,2 (23) 

.1 + {2k) ] 

with A = 2T,/T, and Ti = 157807 K the hydrogen ioniza- 
tion threshold temperature. In agreement with the above 
argument, the effective recombination coefficient at level 
2P is approximately 2/3 of the case B recombination co- 
efficient and that is what we use to convert recombination 
rates into Ly-a photon emission rates. Thus we assume 
that the conversion efficiency from recombination to Ly-a 
photons is exactly the same for all simulation cells. This 
is a good assumption since the conversion efficiency has a 
very weak temperature dependence. 

The exact conversion efficiency for each source cell also 
depends on the rate at which collisions redistribute atoms 
between the 2S and 2P state. Collisions with both elec- 
trons and protons are relevant. To get an idea for the cross 
sections involved, for a temperature of IG^K and thermal 
protons C72S^2P ^ 3 X lO-^^cm^ (Osterbrock 1989). For 
thermal protons and electrons the thermally averaged col- 
lisional cross sections for the processes 

H{2P)+p^ H{2S)+p (24) 

and 

H{2P) + H{2S) + e (25) 

arc Qp — 4.74 x lO^^cm-'/s and qe ~ 5.70 x 10~'''cm'''/s, 
respectively, for a temperature of lO^K (cf. table 4.10 of 
Osterbrock 1989). The 2P to 2S transition is relatively 
important when the proton number densities are small (< 
10*cm~^), and in this case there is some probability that 
the Ly-a photon gets destroyed through a two quantum 
decay. For higher densities the opposite conversion (2S 
to 2P) becomes important, canceling out the destruction 
effect (Osterbrock 1989). At the lower density regime, 
which is applicable in the simulations since there rip < 
10^ cm^^ everywhere (within the virial extent the proton 
number density range is lO"**— 10^'^ cm~'^, with most cells 
in the range 10~^ — 1 cm~^), we can check how important 
this process really is by comparing the radiative decay time 
and the typical time between collisions, 

^ qp{T)np + qeiT)ne 

^ A21 

~ 8.5 X IQ-^^ripT^^-^'' (26) 

where the number densities of protons and electrons were 
assumed to be roughly the same and in cm~^, A21 = 
6.25 X 10^s~^ is the spontaneous radiative decay for the 
Ly-a transition, and temperature is measured in lO^K 
units. The temperature dependence of the collision rates is 
taken from Neufeld (1990). For the temperature and pro- 
ton/electron density ranges relevant to the source cell con- 
ditions, the probability for a collisional 2P to 25 transition 
is negligible, at least for the initial emissivity. We discuss 
their effect during scattering of the photons in §3.5.1. 

One assumption that we make is that the cascading 
of the Lyman series photons, as well as the re-emission 



16 



TASITSIOMI 




-5 5 10 4S 50 

Fig. 8. — Left panel: Cumulative probability distribution of ccntcr-of-liiic optical depths of ART simulation cells within the virial extent. 
The three different lines correspond to the cell optical depth distribution in Ly-a (solid), Ly-5 (dotted) and Ly-limit (dashed) photons. Right 
panel: Ly-a center-of-line optical depth of the simulation cells within the virial extent of the emitter plotted against the cell recombination 
rate. Since only cells with the highest recombination rates (> 10'"'"'" s^^ or, cquivalcntly, luminosities roughly > -5 X 10'''" photons s^-'-) need to 
be used as source cells, and almost all of these cells have r > 10"^, roughly speaking our 'on-the-spot' approximation is satisfactory (see text 
for details). 



and re-absorption of photons from recombination to the 
ground state, is done 'on-the-spot', namely, locally. In 
our case "locally" means within the same simulation cell. 
This assumption is essential if one wants the Ly-a emis- 
sivity of a cell to depend on its own recombination rate 
only. If not, one faces the complicated situation where the 
Ly-a emissivity of one cell depends on the recombination 
rates and photon cascade processes that are happening in 
other cells as well. The validity of our assumption de- 
pends on the optical depth of Lyman series and ionizing 
photons when traversing a typical cell in the simulation 
(and should also be affected somewhat by resolution). In 
the left panel of Figure 8 we show the optical depth prob- 
ability distribution function for Ly-a, Ly-6 and Ly-limit 
radiation. The distribution function has as independent 
variable the optical depth of simulation cells within 10 
physical kpc (~ virial extent) from the center of the emit- 
ter. These distributions are very similar, differing only 
by the values of r because of different oscillator strengths 
and characteristic frequencies. Clearly, in all cases more 
than half potential source cells are not optically thick, and 
this is expected to get worse for ionizing radiation beyond 
the Lyman limit. However, as shown in the right panel 
of Figure 8 the optical depth of a cell correlates with its 
recombination rate. In this figure the optical depth plot- 
ted is that for Ly-a photons, but it is easy to see how 
this scales approximately with optical thickness for other 
Lyman-series photons. Since only cells with recombination 
rates higher than lO^'^ (or cquivalcntly with luminosi- 
ties higher than roughly 5 x 10^" photons s^^) are used 
as source cells, our 'on-the-spot' assumption seems pretty 
satisfactory, if not always accurate. It becomes less and 



less accurate the higher we go in the Lyman series, and of 
course beyond the Lyman limit but for the time being we 
content ourselves with this approximation, given the com- 
plexities introduced when this assumption is not adopted. 
We will investigate this point further in the future. 

Lastly, to get an idea about the physical conditions of 
the highest recombination rate (luminosity) source cells, 
they consist of two classes with respect to temperature 
and neutral hydrogen fraction: one class contains cold gas 
elements (T ~ lO^K), with a neutral hydrogen fraction 
> 0.9 (and high gas number density). The second class of 
very luminous cells consist of warmer gas elements (T ~ 
10^ K and a bit higher). In the context of Ly-a cooling 
radiation, discussed in the next section, the first class of 
cells are unable to cool via atomic hydrogen cooling since 
they are cold, whereas the second class of most luminous 
cells could cool via atomic hydrogen cooling temperature- 
wise, but that is not happening because these cells are 
highly ionized. 

3.2.2. Ly-a photons from collisional excitations 

A collisional emission mechanism whose importance for 
the simulated objects can be assessed relatively easily is 
that of atomic hydrogen cooling. Using the expression by 
Hui & Gnedin (1997) for the hydrogen cooling rate (used 
in the ART simulations analyzed here), and assuming for 
the moment that this energy is all emitted in the form 
of Ly-a photons, we obtain for the luminosity (number of 
Ly-a photons/s) emitted by a cell 

g-1.18355/T5 

LeooJ = 4.6 X 10-8— — — o^neUH/I/ (27) 
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Fig. 9. -- Maximum cooling Ly-a luminosity, Lcoo!; plotted 
against recombination Ly-a luminosity, Lrec, for all ART simula- 
tion cells within the virial extent of a Ly-a emitter at z ~ 8. The 
cooling luminosity is the maximum possible Ly-a luminosity from 
cooling because it is derived assuming that all cooling radiation is 
emitted in Ly-a photons. The solid line shows the case where the 
two luminosities are equal. Since, as discussed in the text, only cells 
with luminosities roughly above 5 X 10^" photons s"'^ contribute 
significantly to the luminosity of the emitter, this figure shows that 
recombination is dominant over cooling Ly-a radiation. 

with T5 the temperature in units of 10^ K. This is com- 
pared with the recombination luminosity Lrec (— 0.68r) 
in Figure 9. Taking into account the results of the conver- 
gence test performed to specify what is the minimum cell 
luminosity that needs to be taken into account (~ 5 x 10^° 
s~^), we see that the cells which are relevant are cells where 
recombination processes dominate, as can be seen in Fig- 
ure 9. Namely, similar to previous studies (e.g., Fardal 
et al. 2001) we find that the cooling radiation Ly-a con- 
tribution is subdominant compared to the recombination 
contribution, hence in what follows we focus only on the 
latter. 

3.2.3. Supernovae Remnants (SNR) 

A Ly-a source that yields Ly-o; photons both from re- 
combinations and coUisional excitations is supernova rem- 
nants (SNR). ShuU & Silk (1979) have computed the time- 
averaged Ly-a luminosity of a population of Type II SNR 
using a radiative-shock code. They find that the Ly-a lu- 
minosity of a galaxy due to SNR is Lsnr = SxlO'^^n]f'^E^- 
ergs/s, with uh the ambient density in cm~^, Eq the typi- 
cal supernova energy in units of 10'"'^, and Nsn the number 
of supernova per year. Strictly speaking, this quantity also 
depends on the assumptions on the IMF, and the lower 
and upper stellar masses of the mass range over which the 
IMF is to be integrated. This expression includes both 
contributions, from recombination and collisional emission 
mechanisms: from UV and X-ray ionization (coming from 
the hot SNR interior) of the surrounding medium and from 



cooling shells, respectively. A thorough investigation of 
the relative importance of SNR Ly-a emission with respect 
to that from young stars photoionization has been carried 
out by Chariot & Fall (1991, 1993). The general conclu- 
sion reached is that for a broad range of physical conditions 
and assumptions, the SNR contribution is at best a factor 
of 2.5 less than that from stellar ionizing radiation. These 
results make the effort to include the (anyway not resolved 
in ART simulations) SNR contributions superfluous. 

3.3. The Ly-a emitter before RT 

To get an idea of the size of the emitting region, the 
prevailing physical conditions, and for comparison with 
results obtained later after including RT, in this subsec- 
tion we briefly present the emission spectrum and image of 
the emitter as they would appear to an observer at z = 
if the Ly-a photons escaped without any scattering. An 
image and a spectrum of the emitter along a certain di- 
rection of observation is shown in the left and right panel, 
respectively, of Figure 10. 

The image is a surface brightness map (in units of ergs 
s~^ cm~^ arcsec"^) of a roughly 1.4 x 1.4 arcsecs^ field 
which corresponds to approximately one third of the virial 
extent of the dark matter halo the emitter lives in (with 
the virial extent ~ 20 physical kpc in diameter). There are 
two distinct emission regions, each one corresponding to 
the two progenitors that merged and formed this object. 
The color scale for the surface brightness is logarithmic. 
Clearly, the emission region is very small (the largest of the 
two structures is at most ~ 2 — 2.5 physical kpc in diam- 
eter, if one includes the faintest pixels), compared for ex- 
ample to the virial extent of the dark halo. The resolution 
of this image is 0.01 arcsecs (~ 0.05 physical kpc), at least 
10 times higher than the best resolution currently avail- 
able. As discussed before, for these results only cells with 
recombination rates higher or equal to 10^^ s~^ are used. 
The initial frequency is chosen according to a Voigt profile 
that is sampled for each cell out to 10 Doppler widths and 
shifted around the bulk (peculiar + Hubble) velocity (;oni- 
ponent along the direction of observation. The number of 
photons used (3 x 10^) has been determined after a con- 
vergence study. Note that when we study the convergence 
with respect to the number of photons we take into account 
that this must be done in parallel with how far away in 
the wings we go when sampling the emission Voigt profile 
of each cell, since the higher the number of photons used 
the better one can sample frequencies further away from 
resonance. The convergence procedure gave the aforemen- 
tioned number of photons and initial emission frequency 
range (i.e., 10 thermal Doppler widths). 

In the right panel of Figure 10, the frequency resolu- 
tion is A/AA ^ 50000. The line shape has converged, 
iVgiJifimely the peaks shown correspond to real velocity sub- 
structure. For example, the most pronounced peak at 
A = 10952A corresponds to the component of the pecu- 
liar velocity along the direction of observation of the most 
luminous pixel of the image shown at the left panel (with 
coordinates on the image (0.24,-0.42) arcsecs, roughly). 
The dominant contribution to this pixel comes from the 
highest recombination cell of the emitter with a recombi- 
nation rate equal to ~ 1.3 x 10^^ s~^ and a peculiar ve- 
locity component along the direction of observation equal 
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to 0.27 X 10^"^ the speed of light. As mentioned, for each 
emission cell the Voigt profile was used and sampled up 
to 10 thermal Doppler widths. The total width of the line 
however is dominated by the bulk velocity structure of the 
emitter. The full width of the line at the minimum flux 
level shown in the figure (10~^^ ergs cm~^ A~^) is 
roughly 15 A (with the width if bulk velocities are set to 
zero being less than half this). This width corresponds 
to projected velocities along the direction of observation 
roughly in the [—200, 200] km/s range (this is just approx- 
imate, note however that the line is not symmetric around 
the rest frame resonance). This velocity range is what 
is expected given the peculiar velocities of the emitting 
cells. Also shown is the spectrum of the smallest of the 
two substructures {dotted line) of the image shown in the 
left panel. One can easily infer what the spectrum of the 
large Ly-a substructure looks like. 

The results discussed in this section may be specific to 
the emitter at hand, but the considerations themselves are 
pretty general. The same kind of procedure must be re- 
peated for each individual emitter identified in the simu- 
lations. 

3.4. The Ly-a emitter- after RT 

It is interesting to first treat the emitter as a finite con- 
figuration. In this case, as soon as the photons exit this 
configuration (whose size is taken to be roughly equal to 
the virial extent of the object, namely 10 physical kpc) 
they travel towards the observer. In other words at first 
we ignore the effect of the GP absorption. This context is 
pretty similar to that of §2.2 and §2.2.3. We focus on the 
emergent spectrum shown with the solid line in the left 
panel of Figure 11 . 

The spectrum converges if 3 x 10'"' photons are used, 
namely if the same number of photons are used as the 
number of photons needed for the initial emission results 
(discussed in §3.3) to converge. Of course, the higher the 
number of photons the better one samples low intensity 
wavelengths. We find that the number of photons used 
affects wavelength ranges with flux less than about 10~^^ 
ergs cm~^ A~^. The spectrum shown in Figure 11 is 
produced using 10^ photons. The spectral resolution used 
in the figure is A/AA ~ 5000, whereas the spectrum is 
identical if ten times better resolution is used. We have 
performed a large set of convergence tests among which 
the most interesting are for different (smaller) cores for the 
acceleration scheme discussed in §2.3.2, and /or a larger 
minimum To<p{xi) for which the acceleration scheme dis- 
cussed in §2.3.1 is used. Our results are pretty robust, 
as should following the discussion in §2.3.1 and 2.3.2 with 
respect to the one cell convergence results. 

Even though meant for a slab, it is interesting to check 
if some predictions of the Neufeld solution, such as the fre- 
quency where the spectrum has a maximum (~ 0.9(q;to)-^/^), 
are roughly in agreement with the spectrum of the simu- 
lated emitter. Of course, the Ly-a emitter environment 
is neither isothermal nor homogeneous, and it is not ob- 
vious how to define an 'effective' temperature and opti- 
cal depth for these purposes. Thus, focusing on order of 
magnitude checks, setting the expression for the frequency 
where the peak emission occurs in a slab equal to the fre- 
quency where the spectrum of the emitter peaks (say in 



red wavelengths) one finds that the 'effective' optical depth 
and 'effective' temperature of the equivalent slab (i.e., the 
slab that would give a spectrum with peak at the frequen- 
cies where the emitter spectrum peaks) roughly satisfy the 
relation t^T ~ 1.4 x 10^" with T measured in eV. The ef- 
fective optical depth will be at least equal to the most 
optically thick cell the photon found itself in. Since the 
emission originates from the most optically thick cells (see 
Figure 8), tq will be at least 10^. If we assume for ex- 
ample a temperatm'c T = lO'^K, the above relation yields 
To — 10^ which is roughly the optical depth from the center 
of the object to its virial radius along the direction of ob- 
servation. Thus, the maximum of the spectrum is roughly 
where it is expected to be if one assumes the scaling from 
the Neufeld solution (~ 2550 km/s). 

The emerging spectrum looks pretty similar to the spec- 
trum that would emerge from a static configuration, namely 
it has two quite similar peaks, one to the red and one to 
the blue of the Ly-a resonance. Note however that the 
peaks are not really symmetric, since the flux decreases 
more rapidly near the resonance. The width of the blue 
peak at a flux level of 10~^^ ergs s~^ cm~^ is roughly 
180 j4 or ~ 5000 km/s. We obtain quite a similar spectrum 
if we set the bulk velocity field to zero in the code, that is 
kinematics do not seem to play a crucial role in this case. 
In the case of the specific Ly-a emitter and for the specific 
direction of observation, analyzing the bulk velocity field 
(i.e., the peculiar velocity field since the Hubble expan- 
sion is negligible at the distances we are working) we find 
that there is some net infalling motion, but with signifi- 
cant transverse velocity components as well. Hence, the 
obtained static like spectrum does not come as a surprise. 
Furthermore, the peak asymmetry due to the existence of 
bulk fields depends on the relative magnitudes of the bulk 
and thermal velocities (e.g., if the bulk velocity is close to 
the thermal we do not expect a significant asymmetry since 
one scattering can give, e.g., a red photon moving in a con- 
tracting medium a large enough shift to erase the effect of 
the contraction) which varies from cell to cell, and it also 
depends on the optical thickness. Since thermal velocities 
are typically small compared to bulk velocities in simula- 
tions, the optical thickness is a more crucial factor. For 
such extremely optically thick media where the spectrum 
is expected to have in the context of the Neufeld solution 
a typical frequency of ~ 2550 km/s, bulk velocities of at 
most some hundreds km / s will not really favor blue versus 
red photons (even if the bulk motion was purely inwards) 
that much, since both red and blue photons see a very 
optically thick medium. 

It would be interesting to have a sense of what is the 
number of scatterings each photon undergoes before exit- 
ing. With the acceleration methods that we have to use 
though it is difficult to keep track of this quantity. A 
simple way to obtain an order-of-magnitude idea of the 
number of scatterings in such or, at least, similar configu- 
rations can be obtained by one of the examples discussed 
in §2.2.3. For the most optically thick case (tq ~ 8.3 x 10^) 
and a point source emitting photons that propagate in a 
stationary medium the number of scatterings in one run of 
~ 2000 photons varies from 2.5 x 10^ up to 4 x 10'', with an 
average of 8.3 x 10®, and a median of 6.6 x 10®. Two thirds 
of the photons are in the [4.6 x 10®, 2.1 x 10^] scatterings 



Ly-a radiative transfer in cosmological simulations 



19 



-17,5 -16,6 -1.5,7 -14,8 -14,0 



U 
0) 

O 



O0.5 - 



- 



-0.5 



4 



0.5 
arcsecs 



o 



17 



°< -U 



19 



-20 



-21 



-22 




1.0945 



.095 1.0955 
A(10'*A) 



1.0965 



Fig. 10. — Left panel: Image of Ly-a direct emission (i.e., assuming the Ly-a photons escape directly to the observer after they are 
produced). The approximately 1.4 X 1.4 arcsecs^ (~ 6.5 X 6.5 physical kpc) field corresponds to roughly one third of the virial extent of the 
dark matter halo where the emitter lives. The surface brightness (SB) is bolometric and in units of ergs s~^cm~'^ arcsecs"^. The SB color 
scale is logarithmic. The object had undergone a recent merger, that is why there are two distinct luminous blobs that dominate the emission. 
Right panel: Initial Ly-a injection spectrum. Shown are the total spectrum {solid line), namely the spectrum for the image shown in the left 
panel, and the spectrum of the smallest of the two blobs in the image [dotted line). Note that the wavelength is in 10*A (i.e., fim). The 
dashed line shows the Ly-a resonance for z ~ 8. See text for discussion of the structure of the line. 



range. More generally, we find that similar to the Neufeld 
problem, the average number of scatterings in this spheri- 
cal configuration scales linearly with optical depth at such 
thick media (see discussion in §2.2.1), with the proportion- 
ality constant of order unity. From this linear scaling of 
the average number of scatterings with optical depth, one 
can obtain a rough idea of the average number of scat- 
terings of photons in the simulation environments (for the 
cell optical depth range in the simulations see, e.g., the 
left panel of Figure 8). These numbers also make clear 
why it is absolutely not feasible to perform Ly-a RT in 
the much thicker and more complicated simulation envi- 
ronments without some acceleration schemes. 

Photons at very optically thick regions have to shift off 
resonance significantly to escape, and hence are the ones 
responsible for the significant line width of the spectrum 
(along with the 1/x^ behavior of the wing optical depth, 
as discussed previously). It is meaningful to ask whether 
one should really care about these photons, or instead ig- 
nore them because may be they are trapped indefinitely 
(for any practical purpose) in the dense cells and do not 
participate in the radiation propagation. To answer this 
question we estimate the photon diffusion time and com- 
pare it to the sound crossing and dynamical time scales 
(other time scales, such as the Hubble time scale for ex- 
ample which is ~ 1 Gyr at z = 8 are clearly large enough 
to be non-relevant). Same as with the number of scat- 
terings, to find the exact diffusion times one should fol- 
low the detailed RT. Given our acceleration methods this 
is not done. Instead we use some useful scalings. Since 
the average number of scatterings in very optically thick 



media is roughly equal to tq, then the diffusion time is 
roughly td — NsJmfp/c with Imfp the mean free path be- 
tween scatterings defined through (t) ~ Te~'^dT = 1. 
In other words, since tq = na{x — 0)L, r — na{x)l, then 
the mean free path in units of the total (half) width of the 
slab is a{x = 0)/it{x)1/to, with a{x = 0) the cross section 
at the line center and (t{x) the cross section calculated at 
an effective x so that the above definition for the mean 
free path is valid. Substituting in the expression for td we 
obtain td ~ (t{x = 0)/a{x)L/c. For a slab with tq = 10^ 
we obtain a mean number of scatterings equal to 9.5 x 10^ 
and a median equal to 7.2 x 10^, whereas 67% of the pho- 
tons have between 5.1 x 10^ and 2.3 x 10^ scatterings. For 
the mean free path we find a mean equal to 2.4 x 10~^, 
a median equal to 1.9 x 10~® and 67% of the scatterings 
correspond to mean free paths between 1.4 x 10~^ and 
7.7 X lO"*", all in units of the (half) width of the slab 
L. For the total distance traveled by the photons before 
escaping, we find an average distance of 40.2, a median 
of 32.3 - implying a (t{x — 0)/a{x) ratio of order 10 - 
whereas 67% of photons exit after traveling a distance be- 
tween 16.7 and 96.3, with these numbers as before in units 
of the width of the slab. Based on spatial random walk 
arguments one would have Ngc ~ Tq, hence the distance 
before escape would be ~ tq or 10® for the specific exam- 
ple we use here. However, as discussed in §2.2.1 Ngc scales 
linearly with tq and this makes a big difference. We find 
that the sound crossing time is significantly higher than 
the dynamical time for most simulation cells, hence the 
latter is the relevant time against which the diffusion time 
must be compared. We find that the dynamical time scale 
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Fig. 11. — Left panel: Emerging Ly-a emission spectrum before adding the Ly-a Gunn-Peterson absorption (GPA) (no GPA, solid line), 
with the GPA but without the red damping wing (GP, no DW, short-dashed line), and with GPA and the damping wing [GPA+DW, dotted 
line). Note that the wavelength is in 10* A (i.e., fj,m). The long-dashed line shows the Ly-a resonance for 2 ~ 8. Right panel: Image of the 
Ly-a emitter after RT, GPA and DW. The 1.4 X 1.4 arcsecs^ (~ 6.5 X 6.5 physical kpc) field corresponds to roughly one third of the virial 
extent of the dark matter halo where the emitter lives. The surface brightness (SB) is bolometric and in units of ergs s~^cm~^arcsecs~^. 
The SB color scale is logarithmic. 



is at least three orders of magnitude or more larger than 
L/c which is within a factor of order 10^ - for the various 
physical conditions in the simulation cells - representa- 
tive of td- Note that this comparison also justifies the use 
of 'static' simulation outputs where the RT is performed, 
even though we plan on investigating the possibility of in- 
corporating the RT scheme into the dynamical evolution 
in the simulations. Furthermore, the effect of the simula- 
tion resolution on these conclusions will be investigated in 
a future study. 

The scattering process diffuses the initial number of 
emitted photons on a larger area and hence lowers the 
nMm6er surface brightness (i.e., number per s cm^ arcsec^ 
rather than energy per s cm^ arcsec^). In general the sur- 
face brightness itself can go either up or down, depend- 
ing for example on the velocity structure of the medium. 
To quantify this effect on a photon-by-photon basis we 
choose to calculate the distance on the plane of the image 
between the initial emission point and the point (pixel) 
where the photon makes its maximum contribution to the 
image (see §2.4 on how spectra and images are obtained). 
We find that these distances vary from roughly 10"'^ to 
10 physical kpc, with a median of 0.27 kpc and a mean of 
0.31 kpc. Given that the largest of the two emission re- 
gions has a diameter of ~ 2 — 2.5 physical kpc (see. Figure 
10), this means that the 'size' of the luminous part of the 
object increases on average by more that 10% due to scat- 
tering. If, instead, we focus on the region where a certain 
fraction of photons originates from we obtain quite similar 
results. For example, ignoring the effects of RT, 90% of 
the emitted photons that would reach the observer would 
originate within a radius of roughly 2.5 physical kpc. The 



same percentage of photons after taking into account RT 
would come from a radius of roughly 2.9 physical kpc.-'^'^ 

So far we have been ignoring the GP absorption. When 
adding this absorption we consider two distinct cases. In 
the first case we include the red damping wing of the GP 
absorption, and in the second case we set it equal to zero. 
The latter best-case scenario is what would happen if for 
example the emitter was inside the HII region of a very 
bright quasar. The spectrum obtained in the first case is 
shown with the dotted line in the left panel of Figure 11, 
whereas the spectrum in the second case is shown with 
the short dashed line. An image of the emitter as would 
appear on earth with the GP absorption and the damping 
wing included is shown in the right panel of Figure 11. 

Not surprisingly, when the damping wing is not taken 
into account the spectrum is identical with that before the 
GP absorption with the difference that all fiux blueward of 
the Ly-a resonance is missing. When including the damp- 
ing wing the maximum fiux is suppressed by roughly a 
factor of 61.7 with respect to the maximum flux without 
it. This line is still quite wide, with a width of approxi- 
mately 1370 km/s at a flux level of 10"^^ ergs s~^ cm~^ 
A-\ and a FWHM roughly 620 km/s. 

Lastly, these results have converged with respect to both 
the number of photons and the radius where the detailed 
RT stops (and beyond which the GP absorption is added). 
More specifically, we find that the number of photons re- 
quired for the initial (no RT) emission to converge (3 x 10^) 

Note that the pixels in the right panel of Figure 11 that give the 
impression of a diffusion of the photons due to scattering possibly 
larger than our ~ 10% estimate, correspond to pixels with practically 
zero number of photons. 
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is enough for the with RT and GP absorption resuhs. And, 
the results also converge if a 10 physical kpc radius is used 
for the detailed RT and beyond that the GP absorption is 
added. Convergence has been checked also with respect to 
the minimum cell initial luminosity considered. Wc find 
that the results converge if the minimum luminosity dis- 
cussed in §3.2 in the context of initial emission convergence 
is used. 

3.5. Some additional physics considerations 

Here we discuss the importance of collisions while the 
photons are propagating, as well as the possible role of 
dust (currently not taken into account). 

3.5.1. Collisions 

While the photons are undergoing scattering, collisions 
should be considered in the following three contexts; (i) 
coUisional redistribution within the n = 2 state; if for ex- 
ample a collision makes the atom go from the 2P3/2 to 
the 2S1/2, then the Ly-a photon is destroyed through a 2 
photon decay of the 2Si/2 state. If instead the collision 
takes it to the 2P1/2, the scattering phase function will be 
different and hence it is relevant in cither case to sec how 
probable the coUisional redistribution is (ii) coUisional de- 
excitation of the n = 2, in which case the photon is lost 
(iii) coUisional broadening of the line, which could cause 
non-coherence in the rest frame of the atom. The RT code 
can take all these processes into account, but here we de- 
velop some intuition as to their importance. In fact, since, 
as will be shown, these processes are in practice negligi- 
ble, the corresponding calculations in the RT code were 
switched off when producing the results presented in this 
study. 

Referring to cases (i) and (ii) , the largest coUisional cross 
sections are for momentum changing transitions (AL = 
±1; e.g., Osterbrock 1989). As discussed already, both 
collisions with electrons and protons are relevant, but pro- 
tons are more significant in case (i), whereas electrons are 
more significant in case (ii). We have already calculated 
the probability per scattering that the 2P — > 25 transition 
of case (i) happens (see equation (26)). The maximum 
value of this probability for the conditions of the simu- 
lation cells is roughly 10~^° (assuming T4 = l,np = 10^ 
cm~^, with the latter being of the order of the maximum 
proton number density of cells in simulations. The temper- 
ature dependence is so weak that it docs not really matter 
what temperature one assumes, for order of magnitude 
estimates). So, unless a photon undergoes 10^° scatter- 
ings, collisions of the type (i) should not matter. The cells 
that are relevant for this are optically thick cells where the 
photons scatter repeatedly. Since as we saw Nsc — tq and 
none of the simulation cells has tq larger than a few times 
10^, collisions should not have a significant impact. Note 
that for most cells the number of scatterings for which col- 
lisions may start to matter is orders of magnitude higher 
than 10"'^" (i.e., what is described above is the worst case 
scenario as far as the effect of collisions is concerned since 
it assumes the maximum proton number density, present 
in very few cells). If these collisions do not matter then 
collisions of type (ii), which have smaller cross sections, 
should not matter either. 

In case (iii), if the atom suffers collisions with other 
particles while it is emitting, the phase of the emitted 



radiation can be altered suddenly. If the phase changes 
completely randomly at the collision times, then informa- 
tion about the emitting radiation is lost and coherence is 
destroyed. In this case, in the rest frame of the atom, the 
line profile is Lorcntzian but the total width is the natural 
width plus the frequency of collisions the atom experiences 
on average. Since the importance of this effect as well is 
determined by a comparison of the radiative decay time 
and the time between collisions (i.e., equation 26), from 
the above discussion it becomes clear that it is also negli- 
gible. 

3.5.2. Dust 

Dust absorbs Ly-a photons. Thus, one would assume 
that dust in the presence of scattering that traps photons, 
coiild have a significant effect, and that this may be true 
even if it is present in small amounts, as is expected to be 
the case for the 2; ~ 8 emitter we discuss (with a metal- 
licity roughly equal to 0.1 the solar metallicity) . Indeed, 
Chariot & Fall (1991) found that only a tiny fraction of Ly- 
a photons escape from a static, neutral ISM even if there 
is a tiny amount of dust present. To include the effect of 
dust absorption in simulations we will have to implement 
a recipe to estimate the amount of dust. Even though 
one can come up with an observationally motivated recipe 
(albeit with unknown applicability at redshifts as high as 
8), we postpone such a treatment for a future study, since 
the main focus of the current study is the Ly-a RT scheme 
(which nevertheless includes the probability per scattering 
that the photon will be absorbed, but this probability is 
currently set to 0). 

However, the Ly-a emitter results we present in this 
study should not be taken as unrealistic, since it is not 
obvious how these results will change if wc include the ef- 
fects of dust. More specifically, many starforming galaxies 
are observed to have significant Ly-a hmiinosities (e.g., 
Kunth et al. 1998; Pettini et al. 2000), and this is usually 
attributed to the presence of galactic winds in these sys- 
tems that allow the Ly-a photons to escape after much 
fewer scatterings than in the static medium case. These 
data seem to support the idea that it is the kinematics of 
the gas rather than the dust content that is the dominant 
Ly-a escape regulator. 

Furthermore. Ncufeld (1991) fomid that under suitable 
conditions the effects of dust absorption may actually in- 
crease rather than diminish the observed Ly-a line strength 
relative to radiation that suffers little or no scattering. 
This would happen for example in a multiphase medium 
consisting of dusty clumps of neutral hydrogen embedded 
within a relatively 'transparent' medium. If most of the 
dust lies in cold neutral clouds then Ly-a photons, not be- 
ing able to penetrate those clumps, will not be affected as 
much by the presence of dust (see also Hansen & Oh 2005). 
Although there is no direct observational evidence to sup- 
port this structure for the ISM (i.e., that dust lies preferen- 
tially in cold, neutral hydrogen clumps, even though the 
dumpiness in the distribution of neutral hydrogen itself 
seems to be established observationally (see Hansen & Oh 
2005, and references therein)), such a morphology of the 
dust and atomic hydrogen distribution could help account 
for the lack of strong correlation between dust content - 
inferred from metallicity or submillimeter emission - and 
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hy-a equivalent width. For example, some dust-rich galax- 
ies have substantially higher Ly-a escape fraction than less 
dusty emitters (Kunth et al. 1998, 2003). In addition, Gi- 
avalisco et al. (1996) found that there is no correlation 
between the Ly-a equivalent widths and the slope of the 
UV continuum, which is a measure of the continuum ex- 
tinction and hence of dust content. 

Another reason why it is not obvious how the results 
presented here will change if we take dust into account, 
is that in the current version of the ART code molec;ular 
hydrogen forms only through the catalytic action of elec- 
trons. When molecular hydrogen formation on grains is 
included in the code, some of what is currently taken to 
be neutral atomic hydrogen will transform into molecular 
hydrogen, hence this effect will decrease the optical thick- 
ness of what currently are the thickest cells. 

4. SUMMARY 

We develop a Ly-a RT code applicable to gasdynam- 
ics cosmological simulations. High resolution, along with 
appropriately treated cooling can lead to very optically 
thick environments. Solving the Ly-a RT even for one 
very thick simulation cell takes a long time. Solving it 
for the whole simulation box, or a significant fraction of 
it, takes unrealistic time. Thus, we develop accelerating 
schemes to speed up the RT. We treat the moderately thick 
cells by skipping the numerous core scatterings which are 
not associated with any significant spatial diffusion, and 
go directly to the scattering that takes the photon outside 
of the core. We use depth dependent core definitions, and 
find that quite large core values can be used. For the very 
optically thick cells we motivate our treatment from the 
classical problem of resonant radiation transfer in a semi- 
infinite slab. We find that with some modifications, since 
the simulations have cubic cells rather than slabs, we can 
use the analytical solution derived by Neufeld (1990) for 
the problem of the semi-infinite slab. With these acceler- 
ating methods, along with the parallelization of the code 
we made the problem of Ly-a RT in the complex environ- 
ments of cosmological simulations tractable and solvable. 
Even though our approach assumes a cell structure for 
the simulation outputs, as is inherently the case in AMR 
codes, the Ly-a RT code we discuss is applicable to out- 
puts from aU kinds of cosmological simulation codes. This 
is true since one can always create an effective mesh by in- 
terpolating the values of the various physical parameters. 

We perform a series of tests of the RT code, and then 
we apply it to ART cosmological simulations. We focus 
on the brightest emitter in those simulations at z ~ 8. 
A first interesting result for this emitter pertains to its 
intrinsic emission region and mechanisms. The emission 
region consists of two smaller regions, each corresponding 
to one of the two main progenitors that merged to form 
the emitter at 2; ~ 8. Both regions are pretty small, with 
the larger of the two having a diameter of 2 — 2.5 physical 
kpc. Furthermore, recombination produced Ly-a photons 
is the dominant intrinsic Ly-a emission mechanism, with 
coUisional excitation and SNR produced Ly-a photons be- 
ing subdominant. The intrinsic luminosity of the emitter 
is 4.8 X lO**^ ergs/s, whereas the injection spectrum (i.e., 
initial emission spectrum) shows significant velocity struc- 
ture. 



After performing the Ly-a RT, but before adding the 
GP absorption, the emitter spectrum obtained resembles 
that of a very optically thick static configuration, despite 
the slight trend for inward radial motions. More specif- 
ically, we obtain the usual double horn spectrum. This 
happens because (i) even though there is some net inward 
radial motion, there are still significant tangential peculiar 
velocity components, and (ii) the optical depth is so high 
that velocities of order some hundreds km/s will not fa- 
vor blue versus red photons (i.e., in order to escape, both 
kinds of photons have to shift off resonance much more 
than the shift because of peculiar velocities, thus none of 
the two kinds of photons is favored in particular because 
of the existence of bulk motions). Namely, the velocity 
information is in fact lost because of the extremely high 
optical depth. The width of the two horns is noticeably 
high (~ 5000 km/s), but in agreement with what is ex- 
pected for the high simulation column densities. The size 
of the emitter increases, since the scatterings disperse the 
photons on a larger area. We find that on the plane of the 
emitter image, a photon on average escapes at a distance 
of about 10% of the initial (before RT) emitter size from 
the point it was originally emitted. 

We include the GP absorption in two different ways: 
without and with the red damping wing. In the first case 
the spectrum is identical to that when the GP is not in- 
cluded, with the difference that now we get only the red 
peak (rather than both the red and blue peaks). This case 
would correspond to the situation where the Ly-a emit- 
ter lies within the HII region of a very bright quasar. In 
the second case, where the damping wing is taken into ac- 
coimt, the red peak is also affected. Its maximum flux is 
supprc;ssed compared to when no damping wing is used by 
roughly a factor of 61.7. The resulting line after including 
the wing is still quite broad with a velocity width of about 
1350 km/s at a fiux level of 10~^^ergs s~^ cm^^ A~^, and 
a FWHM of about 620 km/s. The line is quite displaced 
red ward from the Ly-a resonance, and reach a maximum 
monochromatic flux of 10~^°'^ ergs cm~^ A~^. 

Attempting a detailed comparison with existing obser- 
vations, or discussing detection prospects for an object 
such as the simulated emitter is beyond the scope of this 
study. We have studied only one emitter, and this for only 
one direction of observation since our main goal was to use 
it as an application for the Ly-a RT code. Thus, we do not 
have a large enough and representative simulation sam- 
ple yet. Furthermore, currently the highest redshift where 
a Ly-a line has been observed is ~ 6.6 (Kodaira et al. 
2003)^"'^ and it is not known how different the properties 
of higher redshift emitters are from that of lower redshift 
ones. The most recent report at z = 9 is that of Willis & 
Courbin (2005). This study flnds no detections. The sky 
area coverage is possibly a significant factor contributing 
to this no detection result. Instead, we content ourselves 
here with a simple order of magnitude comparison. The 
intrinsic Ly-a luminosity of our emitter is consistent with 
luminosities reported in literature. For example, the high- 
est Ly-a luminosity of the z = 5.7 sample of Hu et al. 

The detection of a z = 10 Ly-a emitting galaxy was recently 
reported by Pello ct al. (2004) following a color selected survey for 
z > 7 galaxies located behind a well studied gravitational lens clus- 
ter, but the exact nature of this source remains contentious (e.g., 
Weatherley et al. 2004) 
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(2004) is roughly 6 x 10^^ ergs/s. Higher luminosities than 
those have been inferred for Ly-a blobs, rather than emit- 
ters. For example, the most luminous blob in the sample 
of Matsuda et al. (2004) has a Ly-a luminosity of 1.1 x 10*^ 
crgs/s. Most observed Ly-a emitters arc unresolved and 
so is expected to be the simulated emitter. Reported sizes 
for the observed objects are in the ^ few kpc range (e.g., 
Hu et al. 2002). Ly-a blobs on the other hand are quite 
more extended, with sizes ~ 100 kpc (Matsuda et al. 2004). 
The widths of the (lower z) observed lines are typically a 
few hundred km/s, whereas the FWHM of the simulated 
line is roughly 620 km/s. As discussed already, the ve- 
locity width of the ART emitter could be affected by the 
very high H column densities which will drop as soon as 
molecular hydrogen formation on dust grains is taken into 
account. In terms of the detectability, if one adopts the 
present day limit of ground based detections of ~ 10~^^ 
ergs s""'^ cm^^ A^^, clearly our simulated emitter would 
be orders of magnitude fainter. If the emitter is embedded 
within the HII region of a bright quasar, in which case 
the red damping wing will be suppressed, the brightness is 
marginally below the sensitivity of current ground based 
instruments. Note, however, that the prospects of detec- 
tion will be much better for JWST which is expected to be 
able to detect ^ 400 times fainter objects than currently 
studied with ground based infrared telescopes. 
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